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CIRCUIT SIMULATION 

This invention was made with Government support under Contract F33615- 
99-C-291 1 awarded by the U.S. Department of Defense. The United States 
Government has certEiin rights in the invention. 
5 Reference to Related Application 

Priority is claimed to co-pending U.S. Provisional Patent Application 
60/261,033, filed January 1 1, 2001. 

Background 

The present invention relates to the design, modeling, simulation, and 

10 emulation of electronic circuitry. More specifically, the present invention relates 
to numerical time-domain simulation of analog or digital electrical circuits using 
mathematical expressions. Present simulation systems suffer from limitations in 
the kinds and topologies of circuits to which they may be applied..^ The complexity 
of systems to be simulated is also limited in current systems by Various 

15 inefficiencies in the simulation modeling process, including state selection and 

state equation building methods. There is thus a need for further contributions and 
improvements to circuit simulation technology. 

Existing techniques for circuit simulation include two approaches. In one, the 
underlying numerical algorithms are based upon a state-space model of the system 

20 being simulated or analyzed. The system to be simulated is specified either in the 
form of a text file using a well-defined syntax, or graphically using boxes or icons 
to represent computational elements such as summers, multipliers, integrators, 
function generators, and/or transfer functions, for example. General purpose 
mathematical packages operate the models by assembling a state-space model of 

25 the system from the user-supplied description. The time-domain response is then 
calculated numerically using established integration algorithms. 

In the other category of solutions, the system is described as an electrical 
circuit, the fundamental branches of which may include resistors, inductors, 
capacitors, voltage sources, and/or current sources, for example. Other circuit 

30 elements such as diodes, transistors, or electromechanical device^ may also be 

defined as a variation or a combination (sub-circuit) of the fundamental branches. 
In this type of simulation system, the model developer describes the'circuit in the 
form of a list of parameters of the fundamental or user-defined circuit elements and 
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the layout of the circuit. Nodal or modified-nodal analysis techniques are then 
employed to simulate operation of the circuit. In that process, the differential 
equations associated with each inductive and capacitive branch are modeled in 
discrete form using difference equations to relate the branch current and voltage at 
5 a given instant of time to the branch current and/or voltage at one or more 

preceding instants of time. The difference equations for the overall circuit are then 
assembled automatically and solved using established methods to simulate the 
time-domain response. 

It is of note that the second category of systems does not force the model 

1 0 developer to derive the state equations or block-diagram representation of the 
circuit. On the other hand, a state-space model of the overall circuit is never 
established or calculated by the program. Consequently, the numerous (and 
generally more efficient) techniques for analysis of state-space mcJidels cannot be 
applied in such systems. f< 

1 5 More recently, a system has been developed that produces a state-space 

realization for the linear part of a system, the non-linear portion and variable 
parameters being represented as separate or external blocks. The overall system is 
then simulated using the appropriate solvers in the host mathematical package. In 
these systems, however, switches are modeled as resistors with very small or very 

20 large values. The resultmg state model, therefore, may be artificiajly stiff and have 
larger dimension than necessary because of the states artificially introduced by the 
resistor-based switch model. In addition, the system does not incorporate mutual, 
non-linear, or variable inductances, capacitances, and resistances in its standard 
library blocks and components. 

25 Some work on an automated state model generator and circuit simulator 

(ASMG) is reported in O. Wasynczuk and S. D. Sudoff, "Automated State Model 
Generation Algorithm for Power Circuits and Systems," IEEE Transactions on 
Power Svstems . Vol. 11, No. 9, November, 1996, pp. 1951-1956. In this work, 
circuits are specified and analyzed using fundamental branches as shown in Fig. 1. 

30 Each fundamental branch includes a switch, resistor rbr, inductor ^.bh voltage 

source ebr, conductance gbr, capacitor Pbr, and current source jbr. tea,ch parameter 
can be fixed or time- varying, and ideal components can be modeled*by setting the 
remaining parameters to zero. Given the parameters for each branch and the list of 
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nodes that the branches connect, the ASMG generates a state-space model of the 
overall circuit. The state-space representation is calculated and solved using 
numerical methods. If a change in switching state occurs, the state model 
generator then recalculates the state-space model and establishes the appropriate 
initial conditions for the new topology. A disadvantage of this system is that it 
cannot be used to simulate circuits that include loops composed of voltage sources, 
capacitors, and/or resistors. This limitation dramatically hindered the ability of the 
ASMG to simulate high-frequency switching transients of power-electronic-based 
systems. 

Subsequent improvements to the ASMG is the fundamental branches 
illustrated in Fig. 2 using the same notation for parameters as was used in Fig. 1. 
Using this form of modeling, more complicated circuit elements could be 
represented, including transistors, diodes, and thyristors. These improvements are 
discussed in J. Jatskevich, "A State Selection Algorithm for the Automated State 
Model Generator," Ph.D. Thesis, Purdue University, 1999. Additional 
improvements are described herein. 

As used herein, a "spanning tree" over a graph (comprising a set of branches 
and a set of nodes to which the branches are connected) is defined to be a subset of 
the set of branches such that at least one branch in the subset is connected to each 
node from the set of nodes, yet no loop can be formed from the branches in the 
subset. 

Also, as used herein, a "topology change event" occurs when the control 
signal for one or more switching elements causes the switching state of that 
element to change. The nodes to which that element is connected will then become 
(or cease to be) galvanically connected to another portion of the circuit. In most 
cases, diis change affects the state equations for the overall circuit. 
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Summary 

It is an object of the present invention to provide an improved system, method, 
and apparatus for simulating electrical and electronic circuits. Another object is to 
provide a system, method, and apparatus to more efficiently simulate the operation 
of a wider variety of electrical and electronic circuits than has previously been 
possible. 

These objects and others are achieved by various forms of the present 
invention. One form of the present invention is a method, including creating one 
or more data structures sufficient to model an electronic circuit as a collection of « 
(at least two) elements. These comprise zero or more LRV elements, zero or more 
CRI elements, and zero or more switching elements. The LRV elements each have 
at least one of (a) a non-zero inductance parameter Ur, (b) a non-zero resistance 
parameter r^r, and (c) a non-zero voltage source parameter Cbr, but'=neither a non- 
zero capacitance parameter, nor a non-zero current source parameter, nor a switch 
parameter. The CRI elements each have at least one of (a) a non-zero capacitance 
parameter Cbr, (b) a non-zero resistance parameter Vbr, or (c) a non-zero current 
source parameter y^,;., but neither a non-zero inductance parameter, nor a non-zero 
voltage source parameter, nor a switch parameter. The switching elements each 
have a switch state and neither a non-zero inductance parameter, a non-zero 
capacitance parameter, a non-zero resistance parameter, a non-zerp voltage source 
parameter, nor a non-zero current source parameter. A first set of state equations is 
automatically generated from the one or more data structures, and, operation of the 
electronic circuit is simulated by application of said first set of state equations. In 
this method, the collection comprises either (1) an LRV element for which at least 
two of L^,^, rbr, or etr are non-zero, or (2) a CRI element for which at least two of 
Cbr, rbr, oxjbr are non-zero. 

In variations of this form of the invention, the simulating step includes 
producing state output data. In some embodiments, some or all of the parameters 
in the first set of state equations change over (simulation) time as a function of the 
state output data. In some embodiments, some or all of the parai^eters change over 
(simulation) time due to a time-varying parameter of at least one eiqment in the 
collection. 


16410-112/154442 


5 


In other variations of this form of the invention, a second set of state equations 
is generated from the one or more data structures upon the occurrence of a first 
topology change event. In some such embodiments, the generating step simply 
involves modifying only the subset of said first set of state equations that depends 
5 on the one or more switching elements that have changed. In other such 

embodiments, each unique vector of switch states represents a topology of the 
overall circuit, and the method also includes (1) storing the first set of state 
equations in a cache; (2) after a second topology change event, determining 
whether a set of state equations in the cache represents the new topology; and 
10 (3a) if the determining step is answered in the affirmative, using the set of state 

equations that represents the new topology to simulate operation of the circuit after 
the second topology change event; or (3b) if the determining step is answered in 
the negative, building a third set of state equations that represents the new 
topology, and using the third set of state equations to simulate operation of the 

15 circuit after the second topology change event. 

In other such embodiments, the method also includes (1) storing said second 
set of state equations in a cache (2) after a third topology change event, deciding 
whether a set of state equations in the cache represents the new topology; and 
(3a) if the deciding step is concluded in the affirmative, using the set of state 

20 equations from the cache that represents the new topology to simulate operation of 
the circuit after the third topology change event; or (3b) if the deciding step is 
concluded in the negative, building a new set of state equations that represents the 
new topology, and using the new set of state equations to simulate operation of the 
circuit after the third topology change event. 

25 Another form of the present invention is a method including creating one or 

more data structures that together store characteristics of some active branches 
■B"""" that make up a graph of nodes and branches that form a circuit, wherein 
B"""" consists of (i) a (possibly empty) set of inductive branches, each having a 
non-zero inductive component but neither a capacitive component nor a variable 

30 switch state; (ii) a (possibly empty) set of zero or more capacitiv^, branches, 

each having a non-zero capacitive component but neither an inductij/e component 
nor a variable switch state; and (iii) a (possibly empty) set B^ of additional 
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branches, each having neither an inductive component nor a capacitive component. 
^acme partitioned into a first branch set * and a second branch set B"^' , 
where the branches in 5,"^""^ form a spanning tree over B""'"^ , giving priority in 
said partitioning to branches not in over branches in B^. B^' is sub- 
partitioned into a third branch set B^„^ and a fourth branch set B^ , where 
^Lk = ^unT n 5^. A fifth branch set B^^ is identified as the union of (i) 5,^^ , 
(ii) B^ n ' , and (iii) those branches in B^ that form a closed graph when 
combined with B^ . B^ is partitioned into a sixth branch set 5,^ and a seventh 
branch set B,^^ , where the branches in B,^^ form a spanning tree over B^^ , giving 
priority in said partitioning to branches in B^ over branches not in,B^. An eighth 
branch set 5,f„ = 5,^ fl 5^ is identified. A set of state variables, is. selected, 
comprising (a) for each branch of 5,^^ . either the inductor current'br inductor flux, 
and (b) for each branch of B^;^^ , either the capacitor voltage or capacitor charge. A 
plurality of states of the circuit are simulated using the set of state variables. 

In a variation on this form of the invention, the partitioning steps each 
comprise an application of a weighted spanning tree algorithm, such as, for some 
positive numbers wl and wc, (a) for the partitionmg of B'^"'^, a minimum spanning 

tree algorithm is used with weight function cadbj) = J ^'^ . 

[ 0 otherwise 

(b) for the partitioning of B^ , a maximum spanning tree algorithm is used with 

weight function a^bj) = h '^^'"^'^^^ ^> ^ . 

[ 0 otherwise 

Another form of the invention is a system, including a processor and a 
computer-readable medium in communication with the processor, where the 
medium contains programming instructions executable by the processor to: 

(1) build state equations for a first topology of an electronic circuit that has at least 
two switching elements, where each switching element has a swifch|ng state; 

(2) solve the state equations at some time to provide a state outputWector, in 
which at least two elements control the switching states of the switching elements; 
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(3) calculate the value of a switching variable as a function of the state output 
vector, wherein the value reflects whether the switching state of at least one of the 
switching elements is changing; and (4) if the value of the switching variable at 
time ti indicates that at least one of the switching elements is changing, determine a 
second topology of the electronic circuit for time tt and obtain state equations for 
the second topology. 

In a variation of this form of the invention, the programming instructions 
include a state equation building module, a solver module for ordinary differential 
equations, and a switching logic module. The building is performed by the state 
equation building module, the solving and calculating are performed by the solver 
module; and the determining is performed by the switching logic module. In some 
embodiments, the obtaining is performed by said switching logic module, while in 
others the obtaining is performed by said state equation building module. 

In still other embodiments, at some time tj, at least two switclting elements are 
each either rising-sensitive or falling-sensitive switches. Rising-sensitive switches 
change switching state if and only if a controlling element of the st^te vector has 
passed from a negative value to a non-negative value. Conversely, falling- 
sensitive switches change switching state if and only if a controlling element of the 
state vector has passed from a positive value to a non-positive value. In these 
embodiments, the function is the arithmetic maximum of (1) a maximum of all 
elements of the state vector that control rising-sensitive switches, and (2) the 
negative of the minimum of all controlling elements of the state vector that control 
falling-sensitive switches. 

A still further form of the invention is a system for simulating electronic 
circuits having a processor and a computer-readable medium in communication 
with said processor, where the medium contains programming instructions 
executable by the processor to read element parameters and node connection 
information from a data stream. The stream includes at least one switch type 
specification selected from the group consisting of: a unidirectional, unlatched 
switch; a bidirectional, unlatched switch; a unidirectional, latchec^ switch; and a 
bidirectional, latched switch. The instructions are further executable by the 
processor automatically to calculate state equations for die circuit glfeen the states 
of switches specified by the at least one switch type specification. 
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Brief Description of the Drawings 

Fig. 1 is a diagram of a generic branch for modeling a circuit using a prior art 
system. 

Fig. 2 is a diagram of two branches used to model circuits in another prior art 
5 system. 

Fig. 3 is a block diagram of the overall process for simulating a circuit 
according to the present invention. 

Figs. 4A-4C are branches used to model circuit components using the present 
invention. 

10 Fig. 5 is a -block diagram of the run-time computational routines in one 

embodiment of th& invention. 

Fig, 6 is a block diagram detailing the inductive link current calculator for use 
with the routines shown in Fig. 5. 

Fig. 7 is a block diagram detailing the capacitive tree voltage calculator for use 
15 with the routines shown in Fig. 5. 

Fig. 8 is a block diagram detailing the resistive network algebraic equation 
calculator for use with the routines shown in Fig. 5. 

Fig. 9 is a block diagram detailing the inductive network state/output equation 
calculator for use with the routines shown in Fig. 5. 
20 Fig. 10 is a block diagram detailing the capacitive network state/output 

equation calculator for use with the routines shown in Fig. 5. 
Fig. 1 1 is a schematic diagram of a switch element. 

Fig. 12 is a block diagram of the interaction between subnetworks as analyzed 
for use with the present invention. 
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Description of the Illustrated Embodiments 
For the purpose of promoting an understanding of the principles of the present 
invention, reference will now be made to the embodiment illustrated in the 
drawings and specific language will be used to describe the same. It will, 
nevertheless, be understood that no limitation of the scope of the invention is 
thereby intended; any alterations and further modifications of the described or 
illustrated embodiments, and any further applications of the principles of the 
invention as illustrated therein are contemplated as would normally occur to one 
skilled in the art to which the invention relates. 

Generally, the system illustrated in Figs. 3-10 and described herein simulates 
Operation of a circuit in the time domain by collecting component parameters 
(variable and/or constant) and the overall circuit topology, then establishing a 
minimal state space for each active topology (as they are encountered), building 
state equations for that state space, and solving those equations fbr; relevant steps in 
time. As discussed herein, the above process can be very much streamlined using 
the additional knowledge and techniques provided in the present invention. 

The overall simulation process will now be discussed with reference to Fig. 3. 
A description 21 of the circuit, including an identification of constant and variable 
parameters, initial conditions, and the like are fed into state model generator 31, 
which provides state model 23 to solver 33. Solver 33 generates the simulation out 
put 25 for data consumers (such as log files, graphical visualization tools, and the 
like) as are known in the art. Solver 33 also provides continuous state information 
27 to switching logic 35, which determines the state of one or more switches in the 
circuit should be changed. The result of this analysis, switching state 29, is 
provided to state model generator 31. If topological changes to the circuit are 
indicated in switching state 29, state model generator 31 updates the state model 
and passes that updated model to solver 33 for continued simulation. 
FINITE ELECTRICAL NETWORKS 

For the purposes of automated modeling, it is assumed that electrical networks 
can be composed of inductive, capacitive, and switch branch mo(|els depicted in 
Figs. 4(a), (b), and (c), respectively. Using this approach, a wide.^^ariety of 
electrical circuits with different topologies can be modeled by apprc^riately setting 
the branch parameters. Such a modeling approach also assumes that only a finite 
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number of branches is allowed for representation of any particular circuit, and that 
all branches have lumped parameters. Electrical networks satisfying this 
assumption can be modeled using a finite-dimensional state variable approach. 
Thereafter, it is possible to derive a finite-dimensional system of ordinary 
differential equations (ODEs) and algebraic equations (AEs) that would portray the 
dynamic behavior of currents and voltages for any such circuit. It is also important 
to note that only "practical" or "reasonable" circuits are considered, which means 
finite energy, finite current, and finite voltage electrical networks. For such 
networks, all energy storage elements such as capacitors and inductors are allowed 
to store only a finite amount of energy at any instant of time and that energy cannot 
change forms instantaneously. Also, no branches are allowed to carry infinite 
current through their components or have infinite voltage between nodes 
throughout the existence of the network including the instant of time that the 
network is switching between topological stages. These restrictibhs ensure a 
proper commutation (or proper transition) of the network from one topology to 
another. Therefore, a class of electrical networks that permits only proper 
commutation and posses a finite dimensional system of differential algebraic 
equations (DAEs) can be defined as a class of finite networks. Such a class of 
finite electrical networks is denoted as N"''^ where a and are the number of 
state variables in the systems of ODEs for the inductive and capacitive networks, 
respectively, and q is the number of branches used to represent the layout of the 
network. Also, a and 3 are often called the network complexity. Thus, hereafter 
only proper and finite electrical networks from the class N"^^ are considered. 
REPRESENTATION OF ELECTRICAL NETWORKS 

This section will define certain terms with respect to electrical networks fi-om 
the class N"*^ . In particular, a network N can be defined by its topological 
configuration, which is best described by the associated graph denoted as G, and a 
set containing branch parameters denoted as P. Also, a particular switch branch 
may be identified as active or inactive depending on its state. Therefore, in order 
to specify N completely, a topological vector s, which would contaii^'information 
regarding whether each branch is active or inactive, should be addecfe A 
topological state vector s would have ones in places corresponding to all currently 
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active branches, and zeros for the remaining branches which were identified as 
inactive for the current topology. Thus a network of a general kind is an object 
from the class N"*^ and is defined as a triple 

N={G,P,s) (2.1) 
5 It is sometimes convenient not having to deal with inactive branches and assume 
that all of the network branches are currently active. Such an assumption would 
significantly simplify notation without a significant loss of generality. Often, a 
network with all of its branches assumed to be active may be referred to 
augmented or having augmented topology. In these cases, the vector s would not 
10 carry any additional information and, therefore, can be omitted from (2. 1). 

For the associated graph G to be a circuit, it is required that any single branch 
of the overall network N must be included in some loop with other branches so as 
to ensure the existence of a closed path for the current flow. Grapiji^ with such a 
circuit property are referred to as being closed. For consistency, Jtireference 

15 direction (from - to +) is assigned to each branch with respect to the two nodes at 
its endpoints. Also, a network N, in the most general case, may be composed of 
several galvanically independent circuits. Therefore, a directed multi-graph 
denoted by Gj can be associated with network N. Such a multi-graph may consist 
of g closed subgraphs. That is 

20 G',= {Gl,Gl...Gl...Gl} (2.2) 

where each k-xh closed subgraph corresponds to its circuit. The subscript d 
indicates that the graph is directed. Even thought all such subgraphs are 
electrically disjointed they may be coupled through mutual inductances and/or 
capacitances between branches. On the other hand, if ^ = 1, the associated network 

25 becomes a single circuit. Sometimes, for the purposes of analysis, it is also 

necessary to consider the associated undirected graph or subgraph, which can be 
denoted at & and G^, and obtained by omitting the information regarding the 
reference direction of each branch. In a closed undirected graph, any branch must 
be a member of some cycle with more than one branch. . ' , 

30 The directed multi-graph Gj consists of a total of q branches fypm the branch 

set , ' 
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(2.3) 


which, in turn, are connected in some arrangement via p nodes from the node set 


Branches and nodes appear in (2.3) and (2.4) in a definite order which is given by 
their respective subscript indices. Using such a representation, it is possible to 
retrieve a particular branch or a node by referring to the respective ordered set, B 
or A^, with an appropriate index-pointer. 

Therefore, the graph is defined in terms of the branch and node sets as 
= (N, B). In general, the network iV has a changing topology. At one instant 
of time, some branches may be inactive and some nodes may be unused. At 
another instant of time, other branches may switch on and off including different 
nodes, thus defining a new topology. Whenever it is needed, the number of 
currently active branches may be denoted as q' and the number of c&rently used 
nodes as p'; whereupon the number of inactive branches and nodes may be 
expressed q - q', and p - p', respectively. Often, when it is clear what the node 
and branch sets are, only their dimensions may be specified instead of actual sets. 
Also, a spanning tree G*^^ can be associated with each k-th closed subgraph G*, 
and a corresponding forest of such trees = (N, By) can be associated with the 
global graph = (N, B). Note that the set of forest branches By is the subset of the 
global branch set, that is B^ <zB. 

There are several ways to represent a graph on a computer. Some methods 
may take advantage of the sparsity of the interconnection matrix and are, therefore 
more efficient than others in terms of the memory required to store a particular 
graph. A less efficient but algebraically convenient method is to represent a graph 
in matrix form. In particular, a node incidence matrix A/ for the multi-graph has p 
rows and q columns (one row for each node and one column for each branch, all 
ordered). Even though this matrix never has full rank, it can be referred to as a full 
node incidence matrix meaning that it is not reduced to the size of its rank. This 
matrix is conveniently formed from positive and negative node incidence matrices 



(2.4) 
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where the entries of A+ are A+(i,j) = 1 if the y-th branch includes the i-\h positive 
node, and zero otherwise. Similarly, the only non-zero entries of A. are A .(/, j)=l 
if they-th branch includes the i-th negative node. For the networks with changing 
topology, this matrix can be updated for each topology such that for each currently 
5 inactive branch bj the corresponding >th column of A/ is replaced with zeros. The 
transpose of A/ is also known as the adjacency matrix of a graph. For large graphs, 
these matrices are quite sparse and, therefore, may be stored using techniques 
optimized for sparse matrices. 

It is understood in the art that each connected subgraph of the network 
10 contributes exactly one zero-row to Arref. Thus, the rank of A/ is always the 
number of active nodes less the number of connected subgraphs, which can be 
written as 


Using the branch models in Fig. 4, the parameters of the circuit can be 
conveniently arranged into parameter matrices denoted as Rbr, Gbr, Lbr, and Cbn In 

20 one embodiment, these square matrices would have dimensions equal to the total 
number of branches in the global network to be modeled, and would contain the 
parameters of each branch in the corresponding diagonal entry of each matrix. If 
some branches are coupled through mutual inductances or mutual capacitances, 
then those mutual parameters are represented in the off-diagonal enti-ies of hbr or 

25 Cbr correspondmg to these branches. In the same way, mutual resistances or 
mutual conductances can be modeled between branches. Even though such 
quantities might not have physical meaning similar to mutual inductances or 
mutual capacitances, they can be employed for simulation purposes. If some 
branch parameters are not present in the circuit, the corresponding entries in the 

30 parameter matrices are filled with zeros. Similar to parameter matrices, the 

external independent voltage and current sources are also represented as vectors 
denoted as ebr and ji,r. These vectors would have voltages and currepts in those 
entries corresponding to the branches with respective sources and zeros elsewhere. 
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In some circuits to be simulated, the network parameters will not depend on 
the currents and voltages applied to the branches. In other circuits the parameters 
will vary with time, so the parameter matrices would also become matrix functions 
of time RfcrCO. GbAf), ^bAf), and CbAf). For simplicity of notation herein, such 
time dependence will not be written explicitly. Also, since inductors and 
capacitors are energy storage components, in addition to their time dependence, 
their total derivatives with respect to time must also be known. Thus, the 
derivatives of time-varying inductances and capacitances are additional inputs into 
the model of the energy transformation and exchange processes in the circuit. 

The network parameter matrices and the vectors of voltage and current sources 
can be grouped to form a parameter set for the global network, which can be 
defined as 

i'={R..L..^.e.,G..c..^,j4:^. a.io) 

It is also convenient to define the following subsets of P 

dL, 
dt 

^c={G,.,C,,^.j,| (2.12) 
P^={R,„G,„e,,,j,J (2.13) 

Clearly, 

P = P^kjPcUP^ (2.14) 
Subdivision of the parameter set P corresponding to the global network N into 
subsets (2.11), (2.12), and (2.13) has the following goal. First, consider a grid of 
branches of either of the two types such that all parameters can be placed in one set 
Pa. a network of this kind does not have energy storage element such as inductors 
or capacitors. Modeling such a network would require solving a system of AEs 
relating branch voltages and currents to external sources through the network 
parameters and topology. From the nature of the equation that needs to be solved 
in order to obtain all branch currents and voltages, such networks ^'referred to as 
being algebraic and denoted as Na- On the other hand, networks whipse parameters 
are placed in the set P/, have inductive elements which may store energy in the 


R...L...-^,e,4 (2.11) 
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magnetic field. Therefore, a network of this type possesses a system of ODEs, or 
more precisely a state equation, whose natural state variable may be inductor 
currents or flux linkages. Similarly, networks whose parameters are placed in Pc 
have capacitors that may store energy in an electric field, and therefore, posses 
5 ODEs whose natural state variables may be capacitor voltages or charges. These 
two networks are referred to as inductive and capacitive, respectively. 

Thus, three types of elementary electrical networks, namely inductive, 
capacitive, and algebraic, denoted as Nl, No andiVA, respectively, are considered. 
There, only the first two networks are allowed to have energy storage circuit 

10 components, and therefore, only these two networks can have state variables. In 
the most general case, an electrical system may be composed of all types of 
branches considered above. Therefore, a corresponding network N can be 
represented as an interconnection of Nl, No and Na, which is symbolically 
expressed in terms of imion as f, 

15 N^Nl^Nc'^Na (2.15) 

Depending on a particular circuit to be modeled, some networks may or may not be 
present in (2.15). 

TOPOLOGICAL CANONICAL FORM 

Models of electrical networks must obey all laws of circuit theory, and in 

20 particular, Kirchhoff's current and voltage laws, KCL and KVL. Automated 

modeling requires formulation of equations describing the network using KCL and 
KVL that are written algorithmically for the entire circuit or its sections. 
Moreover, KCL and KVL may be written for a network of the most general kind 
based only on its topological configuration and regardless of the actual branch 

25 elements of their volt-ampere characteristics. Therefore, the starting point in the 
analysis of graph G| = {N,B) is the associated full node incidence matrix A/. 

It is necessary to express KCL and KVL for the network whose topology is 
given by Gj = {N,B) and stored in A/ For simplicity of notation, all branches in 
this section are assumed to be active. The procedure continues as follows. First, 

30 instead of Gj = {N,b) , the corresponding undirected version G*.== (N, B) is 

considered. Then, suppose that from G* = (AT, B), it is possible to fi*id a forest of 
spanning trees G*„, = {N,By )that spans the entire node set iV using the set By, 
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which is some subset of B of minimal size over all branch subsets. In general, 
^trees = ) ^0"^ being unique for a given = (A^, B). In fact, it can be 

shown that a graph whose full node incidence matrix is A/ has det(A/A/) total 
number of spanning trees. Moreover, the partitioning of a multi-graph into 
spanning trees G*„, = [N,By ) and remaining link-branches is a nondeterministic 
topological problem. One way to simplify this task is to associate appropriate 
weights with the branches and convert the problem into finding a spanning tree 
with the minimized/maximized sum of such weights. This approach is known in 
network optiniization as the minimum/maximum spanning tree problem. There are 
several well-known minimum/maximum spanning tree algorithms which roughly 
yield the order of complexity Oiq log p) and better depending on the data structure 
used to represent a graph. Also, the branch weights should be relatively simple so 
as to promote good performance as well as to be able to prove certain useful 
properties. Thus, the weights can be assigned to network branches' based on their 
respective parameters (2.10). This method of obtaining a spanning tree and a set of 
links, with some desired property, will be utilized in the present, exemplary 
embodiment for the purpose of automated modeling. 

For the purpose of this section, any proper forest of spanning trees would 
suffice. Having performed such graph partitioning, the set of tree-branches By is 
identified. Thereafter, the set of remaining link-branches Bx can be determined as 
B,=B-By (2.16) 

Since the subsets By and Bx are identified based on G*^„ = (N,By), it is now 
also possible to re-sort the complete branch set B such that all tree-branches appear 
first from the left and all the remaining link-branches on the right. That is, the 
branch set B can be re-ordered as 

B = {By,Bj (2.17) 
The new branch order defined in (2. 17) can be applied to the columns of A/. 
In this case, the new branch order can be related to the original branch order 
through a permutation matrix which is denoted as Tp. Multiplyiiig^from the 
right by Tp results in a matrix whose columns are ordered to corres^nd to the 
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branch set B in (2.17). That is, the permutation matrix Tp should sort the columns 
of A/ such that 

A;.T, =[a,_.A;^] = A^ (2.18) 
This permutation matrix can be assembled from an identity matrix by sorting its 
columns at the same time as the branches in (2.17). Note that the multiplication 
from the right by Tjperfonns the reverse column permutation and restores the 
original column order. That is 

A^=A/T; (2.19) 
Here, as well as further on, the hat sign above the matrix denotes that the 
corresponding matrix or vector quantity is referred to a branch order different that 
the original order. 

After the full node incidence matrix is expressed in the form (2.18), its 
algebraic properties can be utilized. In particular, since the colurhins of A,r«e, 
correspond only to branches that form a free-of-cycles forest, this matrix must have 
a full rank. On the other hand, Aunks contains only columns corresponding to the 
link-branches, and therefore, none of its columns add to the rank of A/. Therefore, 
the rank of A/, as well as the rank of A/, is determined by the columns of A.,rees- 
Therefore, these columns can be chosen to be the basis columns of A/ or A/ With 
respect to (2.18), the basis colunms appear first from left-to-right. The RREF of 
(2.18), therefore, will have die identity matrix of the size of its rank as the upper- 
left block, and its re-ordered structure becomes 


I„,A^ 


(2.20) 


The zero rows on the bottom of (2.20) do not include any useful topological or 
algebraic information about the network and therefore can be deleted. We define 
25 the topological canonical form (TCF) of the node incidence matrix to be the result 
of the reduction of (2.20): 

A,,,=rCF(A,) = LA^J ' (2.21) 

The TCF is one of the key concepts used herein to describe netv^^orks. From 
Atcf, the reduced incidence matrix A„ and the so called basic loop^itoatrix are 
30 obtained as 
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^a=^TCF (2.22) 

B.=[-ALI..J (2.23) 
Matrices (2.22) and (2.23) were assembled with respect to the branch order (2.17) 
and can be easily recovered for the original branch order as 

A, =A,T; (2.24) 

B, =B,T; (2.25) 
Denoting vectors of branch cvurents and voltages as ibr and \bn respectively, KCL 
and KVL written for the whole network have the usual form 

A„W = 0 (2.26) 
B/,Vft, = 0 (2.27) 
Also, instead of using recovered matrices (2.24)-(2.25), it is possible to re-sort ibr 
and ybr into i^, and v^^ for the new branch order (2. 17) using the ^june 
permutation matrix Tp, and then use these with (2.22)-(2.23) to form KCL and 
KVL similar to (2.26)-(2.27). That is 

AA.T,=AA.=0 (2.28) 
A,v,,Tp=B,v,,=0 (2.29) 
Furthermore, based on the structure of TCF, an interesting property of 
matrices (2.22) and (2.23) can be utilized. This property relates voltages and 
currents of the branches corresponding to subsets By and Bx in (2.17). In particular, 
based on the branch order (2. 17) it is possible to define the vectors of tree and link 
currents and voltages as 

L=[i,'iJ (2.30) 
n.=[v,.vj (2.31) 

Then, based on (2.22), (2.23), and (2.30), the vector of re-sorted branch currents 

can again be expressed in terms of the link currents as 

K=(^bfK (2.32) 

Similarly, the vector of re-sorted branch voltages can be expressed in terms of the 

tree voltages as ; 

^br=(Ky^y (2.33) 
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The relations (2.28)-(2.33) play an important role in formulating the governing 
DAEs describing electrical networks. 
NETWORKS WITH STATE VARIABLES 

In this section, a class of finite electrical networks N^^'^ is considered. In 
particular, the two types of networks that can be modeled using a state variable 
approach are inductive and capacitive. Equipped with techniques based on a 
topological search for an appropriate forest of spanning trees, the conditions under 
which a corresponding state equation can be assembled are also considered. 
Inductive network state equation formulation 

An inductive network can be built using branches of the type depicted in Fig. 
4(a). Assuming that all switches are active, a network of this type is given as 
(2.80) 

f/" (2.80) 

where the parameter set Pl is defined in (2, 1 1). ' 

For a given topology of A'i^ the branches must be re-sorted into subsets as in 
(2. 17). The subset By would collect branches that form spanning trees for all 
subgraphs. Each such tree is free of cycles, and covers all active nodes in its 
subgraph. The second category, denoted as subset B^, takes the link-branches. 
These branches are the links in a sense that addition of any of them to the spanning 
tree would result in a cycle. These branches can carry state variables- 
independent currents — and since their number is minimal for each spanning tree, 
they form a minimal set of states for the Nl- Thus, in order to partition = (N, B) 
and obtain the required branch sets, the following weight function wdib) is defined 


For.-=,...,. w,(^,={°'.ff-''<^--^^:° (2.8:) 

Applying the minimum spanning tree algorithm (MinSTA) to G^ = (N,B) with 

weight function WL(b) a minimum spanning forest G^^^^ = (N,By) is obtained as 
before. Thereafter, the set of link branches Bx can be determined as 5n (2.16). 

■ > 
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An important observation is that the network Nl has no non-inductive loops if 
and only if there is a set Bx in which the number of branches equals the sum of 
their weights 

X = size(B^) = s (2.82) 

Condition (2.82) in necessary and sufficient. The set need not be unique, but 
any such set Bx, satisfying (2.82) is equivalent in a topological sense. That is, an 
arbitrary set of branches (that can be larger than J5J may be chosen to represent 
state variables, such as independent currents, in Nl if and only if it contains a set Bx 
satisfying topological condition (2.82). If however, (2.82) cannot be met, it 
follows that the given network is more than just a single inductive network. An 
algorithm for handling more than one network will be presented in later sections. 

Applying the MinSTA with (2.81), assuming (2.82) holds, the order in which 
all branches are grouped and sorted according to (2.17) is obtained. This final 
branch order is related to the original branch order through a pemiutation matrix, 
which in this case is denoted as Tl, Using this permutation matrix and the TCF, 
matrices Aa and B^, are found as usual. Thereafter, a vector of state Variables for 
Nl may be chosen to be a vector of independent currents, such that 

hr = B^i^ (2.83) 

The corresponding state equation is obtained using the dimensionality 
reduction procedure discussed in Wasynczuk and Sudhoff. The procedure may be 
as follows. The voltage equation written for the network is multiplied from the 
left by the corresponding matrix of KVL and all vectors of branch currents are 
replaced with (2.83). The result is 

f =-I.-(R..^»)l,-I.-B,e,. (2.84) 

Capacitive network state equation formulation 

A capacitive network with augmented topology can be defined' using the 
corresponding parameter set and an associated graph as i 

Nc-iG,Pc) (2.85) 
where the set Pc is defined in (2. 12). 
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Also, in the case of Nc, the forest By of the associated G^, = (N,By) can 
represent state variables-the capacitor voltages. Thus, a slightly different 
topological approach might be used. In particular, a MaxSTA with a different 
weight function can be applied in order to find maximum spanning trees. This 
time, the weight function wdbj) is defined such that 


fori = wcibj) - 


0. if C,,O-,i) = 0 ^2.86) 

1. if C6,o;j)^o 


In this case, the network Nc should not have non-capacitive tree-branches. Such a 
condition is satisfied if and only if there is a branch set By for which the following 
is true 

^ (bf ) = sizeiBy ) = r (2.87) 

Condition (2.87) is also necessary and sufficient, and therefore, the discussion 
thereof applies to this case as well. Again, the set By need not be unique, but all 
such sets for which (2.87) holds are equivalent in a topological sense. If the 
condition (2.87) cannot be met, it can be shown that the corresponding network is 
not just Nc but a union of the form Nq u A''a. However, for the purpose of this 
section, a single capacitive network is considered. 

Thus, applying the MaxSTA with weight function wc(b) and condition (2.87), 
a similar permutation matrix Tc, can be formed such that multiplying A/ from the 
right by Tc the desired branch order (2. 17) for the Nc can be obtained. Further, 
matrices Aq, and B^, are formed exactly as before. 

The natural state variables for Nc are the capacitor voltages. Therefore, a 
vector of state variables Vycan be chosen to be a vector of independent capacitor 
voltages such that 

Vbr = Ajyy (2.88) 
The state equation is very similar to (2.84) and is also obtained using the 
dimensionality reduction procedure appUed to the current equation written in 
matrix form for Nc- The result is ' 'y> 
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where notations are analogous to those used in Wasynczuk and Sudhoff, supra. 
Specifically 

% " K^hr^a (2.90) 
= A„Cft,Aj (2.91) 

Transformation of Variables 

An important observation regarding the selection of state variables for as 
well as for A^c is that, once the set of states in (2.84) and (2.92) are obtained based 
on topological conditions specified previously, none of the actual states are 
required to be associated with any particular branch. That is, a tppologically 
proper set of state variables in Nl and A^ccan be transformed into igi^ equivalent set 
by any non-singular coordinate transformation. Such a transforma:^ion of states is 
nothing more than a change of variables in the state space. To show this, it is 
useful to introduce time-varying coordinate transformation defined by two non- 
singular square matrices K/,and Kc of dimensions corresponding to the number of 
states in iVt and iVc, respectively. Multiplying (2.26) and (2.27) from the left by 

and Kj^ respectively, KCL and KVL still hold: 

^cKhr = 0 (2.99) 
^L^b^hT = 0 (2.100) 

Thereafter, the new state variables for Ni, and Nc are related to \x and Vy as 
= ^L^L (2.101) 
= Kpx,. (2.102) 

And die corresponding branch currents in Nl and the branch voltages in iVc are 
expressed as 

iftr = BjK^Xi = Bftx^ (2.103) 

J. -r ' ' ■ 

Vftr = A^KpXj; = AaXc (2.104) :■ 
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Since the change of variables for the two networks is very similar, only the 
corresponding state equations for Nl will be written explicitly. Proceeding as 
follows 

= -Ki'L;'(R, + £'*]KiX^-Ki'^(K^)Xi-Ki'L;'Bje,, (2.105) 

5 Choosing to be a constant matrix, the result of transformation is 

= -K2'L:'(R,+f jK^x^-K£'L;'B,e,. (2.106) 

where the reduced quantities with the subscript "jc" and tilde above are defined 
using from (2.106) instead of just Since (2.84) and (2.106) are related 
through the similarity transformation via constant matrix Ki, the eigenvalues of the 
10 resulting dynamic matrices in both state equations are the same. Also, since it is 
possible to choose K/, to be a permutation matrix, it follows that the branch 
ordering does not change the eigenvalues of the system. 

It is also interesting to observe what happens to (2.105) if the time- varying 
coordinate transformation is chosen to be 

\ Ki = L;' (2-107) 

In this case, the new state variables become fluxes 

>'x = K24,= LA (2.108) 
Then, with the following equality 


20 equation (2. 105) reduces to state equation 


(2.109) 


_* = -R^-\-Bie2,, (2.110) 

A similar change of variables can be performed on Nc. The equations can be 
simplified even further by considering only constant capacitances, wherein a time 
invariant coordinate transformation Kc can be used. Therefore, (2. 105) may be 
rewritten as . 
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= -Kc'c;'g,Kcx^+Kc'c;'aj,, (2.111) 

= -Cy^Gy+CyAaJbr 

where, as before, the reduced quantities with the subscript "y" and tilde above are 
defined using from (2.111) instead of Aq. 

Even though Kc in (2.100) can be any non-singular matrix, it is always 
5 possible to choose it in a way that the new state variables have particular physical 
significance. Siniilar to Ni^ where the states can be transformed from currents i;cto 
fluxes Xx, in the case of Nc the states may be transformed from voltages v^, to 
charges q^by an appropriate choice of transformation. In particular, by defining 
the transformation as 

Kc = C;' (2.112) 

10 ""y ^ ^^'"y (2.113) ■ 

The state equations (2. 1 1 1) with capacitor charges being states becomes 

J^ = -GyC;'qy+A„jft, (2.114) 

MULTIPLE NETWORKS 

An algorithm for automatically deriving a system (or systems) of DAEs for 

15 time-domain simulation of many practically useful electrical networks will now be 
discussed. As noted previously, there are some conditions under which a global 
network with state variables, here denoted as 

N = (G,P^), where C - C (3.1) 
cannot be structurally represented entirely only as Nl or as Nc. Specifically, an 

20 attempt to represent the network N structurally entirely as Nl fails if there is a 
shortage of inductive link branches in the set B^, On the other hand, the network N 
cannot be viewed only as Nc if the corresponding forest of spanning trees lacks 
capacitive branches in the set By. In these cases, the global network cannot be 
viewed entirely as a single type of a network, and therefore, a more general 

25 approach is required, i • 

A method of handling networks with arbitrary topologies is presented here 
based on the separation of N into interconnected networks. That is, ,if it is not 
possible to represent an entire circuit as a network of a single kind, it is necessary 
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to identify some cuts of N that can be grouped together to form several networks 
based on their topological properties. Separation of N into Ni^ Nc, and Na in a 
way such that it is possible to obtain consistent state equations for Nl, Nc, and a 
system of algebraic equations for iV^i is, therefore, a generalization of the ASMG 
5 approach for circuits with arbitrary topologies. In this sense, a robust state 
selection algorithm is a key to such a generalization. 

The selection of branches that can carry state variables in Nl or Nc is 
determined through a topological partitioning of the global graph & = (N,B) into 
a forest of trees and remaining link branches as 
10 G* = WB) = {G^.G^} (3.2) 

where 

G^ = (iV.B^) (3.3) , 

= (N,B^) (3.4) ' 

In previous discussions devoted to networks with state variables, the actual 
partitioning of = (N, B) was performed with two distinct objectives. First, with 
15 the help of the MinSTA and the weight function WL(b), it was possible to reorder 
branches such the that set of link-branches Ex could be chosen to represent state 
variables in Nl. Second, using a similar topological approach based on MaxSTA 
and the weight function wc(.b), it was shown that branches from the set of tree- 
branches By can represent state variables in Nc. Thereafter, the conditions for the 

20 state equations to be complete were: for Nl the corresponding set B^ must contain 
only inductive branches; and for Nc the corresponding set By must contam only 
capacitive branches. These two branch sets could be obtained independently 
running MinST and MaxST algorithms with different weight functions. These 
topological conditions are stated in (2.82) and (2.87). 

25 The technique of network identification and partitioning introduced here is 

based on the TCF of the node incidence matrix assembled for the global network. 
As it will be shown, the TCF makes topological information about N available in 
an "algebraic" sense for the further use in KCLs and KVLs. There are also some 
details about notation that are worth pointing out. The topological q'jaantities 

30 referred to a particular branch order, such as in (2. 17), will be distin^ished by the 
hat sign, keeping in mind that it is always possible to transform them'back to the 
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original order using a corresponding permutation matrix. The superscripts "L", 
"C, and "A", and the combinations thereof would be employed to relate variables 
and quantities to the inductive, capacitive, algebraic, and the overlapping networks, 
respectively. 

Starting with Nl, in will be shown how an algebraic network can be identified 
and how the corresponding system of algebraic equations can be assembled. It will 
also be shown that even in the presence of Na, a minimal and consistent state 
equation for Nl can still be obtained following the same dimensionality reduction 
procedure set forth in Wasynczuk and Sudhoff. Then, similar derivations will be 
repeated for Nc, in a somewhat simplified form, using its structural duality with 
respect to the inductive case. Finally, a way of obtaining a consistent system of 
DAEs relating all networks will be presented using the established framework. 
Inductive and Algebraic Network Interconnection 

'',i'i4 

Here, it is assiraied that a network iV = fG, Pl) is constructed tising branches 
shown in Figs. 4A and 4C. This time, it is assumed that the MinSTA with weight 
function wtib) is applied, and that in the end (2.82) does not hold. This implies 
that there is no way the set of links Bx can be chosen to contain only inductive 
branches, and 

X wi,{bj) = kKsizeiBJ (3.5) 

It follows from (2.17) and (3.5), there is at least one non-inductive branch in 
addition of which to G,rees = (N, By ) would result in a cycle. This cycle would be 
entirely composed of non-inductive branches. Otherwise, the spanning tree 
corresponding to this link branch would not be minimal which, in turn, would 
contradict the result of the MinSTA. In general, addition of any branch from Gii„ks 
= (N, Bx.) to Gtrees = (N, By) results in a cycle. Therefore, the global graph = (N, 
B) has as many non-inductive cycles as there are non-inductive branches in the 
branch set B^. This property is based on the MinST and the weight function (2.13). 
Thus, based on (3.2), the total number of non-inductive cycles in G^ = (N,B) can 
be determined as 

m = size(B^)-A = s-A (3,6) 

These cycles do not have state variables such as currents or fluxes, ,1'herefore, 
none of the branches participating in such cycles can be placed in Ni^ 
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There are many methods that can be used to identify cycle(s) of a particular 
kind in a graph. Recall that the branches can be reordered similar to (2.17) and 
(2.18) by a corresponding permutation matrix obtained from the MinSTA. This 
time, it is necessary to reorder branches in (2.17) and columns in (2.18) in an even 
more sophisticated way. Specifically, since it is known which branches in B^c are 
non-inductive, they can be identified and put first on the left side in Aii„ks,. Then, 
the branches in By that correspond to cycles linked by non-inductive branches in 5^ 
are identified, and the block A,„es is resorted such that these columns appear on the 
right. Similar to (2.18) this final branch ordering can be expressed in terms of the 
node incidence matrix and permutation matrix as 

A/Ti = CALs.A^.Atk„ALcs] (3.7) 

According to the system of notation employed, the superscript "L" denotes 
that the corresponding branches can be safely placed into inductiv^^network A^£, 
and the "A" identifies all other non-inductive branches, as viewed from Nl, that 
belong to Na. Taking the RREF of (3.7) and removing zero rows on the bottom 
yields a TCF similar to the one introduced in (2.2 1). 


Aj-CT - 


AnxA 


0 I„x„ AJixot! CjixA 


(3.8) 


where 

h- is the number of inductive link-branches in the set Bx, as defined in (3.5); 

m - is the number of non-inductive link-branches in J5» as defined in (3.6); 

T] - is the number of tree-branches in the set By that are linked by the h 
inductive link-branches in B/, 

|j, - is the number of tree-branches in By that are linked by the m non-inductive 
link-branches from Bx. 

Also, Ti + = r, and m + h = s, where s and r are defined in (2.8) and (2.9), 
respectively. The corresponding subsets of tree-branches in G,rees =. (N, By) can be 
readily identified from the node incidence matrix A/ in its TCF such as (2.21). 

The reduced node incidence matrix and the basic loop matrix for the global 
network N are found as usual ' >' 
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(3.9) 


(3.10) 


Even though the TCFs (2.21) and (3.8) are different only in terms of the 
relative order of some columns (branches), expression (3.8) has more advanced 
structural properties that can be utilized in writing KCL and KVL for the sections 
of the global network. Specifically, KCL for A^^ can be written based on (3.9) as 


(3.11) 


Based on (3.10), KVL for AT^ can be derived as 


It is interesting to observe that KCL (3. 1 1) is a self-contained equation for Nl- 
10 Similarly, KVL (3. 12) contains only quantities relevant to Na- In this sense, these 
two equations are de-coupled. 

The vector of state variables for Nl, which is the vector of independent 
currents, is then chosen as 


■[t:r]'--R'' 


In proceeding further, two more equations describing both networks are 
needed, specificaUy KVL for A^£, and KCL for A^^. First, from (3.10), the 
following KVL can be written 

\-{k^.nf {d^i^.f 0,^„ I,,,]v*. = 0 (3.14) 

which may be rewritten using notation for the two networks as j 

~L~L 'LA- A 

Rb-Vhr = Cfc VI, ' 
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and 


Similarly, from (3.9), KCL can be written as 

'A' A "i^'r. 'A, 
Aalbr = Co Ibr = -C 1^ 

5 It is worth noting that (3.15) and (3.17) express the connection between the 

two networks. To be more specific, KVL (3.15) fov Nl now has an additional 
voltage source term that comes from the interconnection of iVi and Na. Similarly, 
KCL (3.17) written for iV has current source term which comes from Nl. 

The dimensionality reduction procedure for obtaining state equations for Nl 
10 remains the same. That is, the voltage equation written in matrix form for the N is 
multiplied from the left by the corresponding KVL loop matrix, ai^tf the vector of 
branch currents is substituted with an expression of the form (2.83). Thus, using 
B/ from the KVL (3.15) and the vector of state currents chosen as in (3.13), the 
state equation for Nl becomes 

where the reduced quantities with the subscript ";c" can be obtained as 

= B^RliBif (3.19) 

L, = BlL.l(Bif (3.20) 


R6r = T^Rj^Ti (3.21) 
Lbr = tJl.^T^ (3.22) 

It can be verified that L;^ is non-singular, and that (3.18) is infaet a minimal 
state equation for A^l with the total of h state variables. Also, note tjjie additional 
source term on the right side of (3.18), which distinguishes (3.18) frdm (2.84). 
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It appears that in order to solve (3.18), the branch voltages for Na should be 
known. These voltages are functions of internal topology, parameters, and external 
sources. In general, it is necessary to compute both currents and voltages for all 
branches in ^Va- Here, KVL (3. 12) and KCL (3.17) are utilized. In particular, 
suppose the currents of Na first, and then the branch voltages. However, in (3.17) 
there are fewer equations than branches in Na- Specifically, (3.17) provides ji 
equations with m + p. unknowns. That is, there should be m more equations, 
precisely one for each cycle in Na. On the other hand, Na has its own voltage 
equation, which may be expressed in terms of the branch order given by the 
columns of (3.7) as 

V6 r = R6rl6r + Cftr (3.23) 

Substituting (3.23) into (3. 11), the following result can be obtained 

Bb Rbrlbr = -BbBbr (3.24) 

where e^^ is the vector of voltage sources corresponding to branches in Na and the 
branch order given by (3.7). 

Combining (3.17) and (3.24), a complete system of /n + ja equations for the 
branch currents of A'^ can be assembled in the following form 


If the branch resistances do not depend on currents, the system of equations (3.25) 
20 has the form Ax = b. In addition, if Na has time-invariant resistive branches, A is a 

constant nonsingular matrix. Therefore, for the time-invariant case, (3.25) can be 

solved by inverting the corresponding matrix A once for each new topology. 

Thereafter, the corresponding vector of branch voltages v^, is found using (3.23). 

Finally, (3.18), (3.23), and (3.25) form a system of DAEs describing the two 
25 networks Nl and A^^. 

Capacitive and Algebraic Network Interconnection 

In this section, it is assumed that a network = (G, Pc) is constructed from 

branch models shown in Fig. 4. It is assumed that the MaxSTA witii the weight 

function wdb) is applied and that condition (2.87) does not hold. This result 
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implies that the set of tree-branches By cannot be chosen to contain only capacitive 
branches, and 


Then, based on (2. 17) and (3.26), the total number of non-capacitive tree- 
5 branches in the forest of maximum spanning trees can be determined as 


When applying the MaxSTA, a permutation matrix that sorts columns of the 
full node incidence matrix A/ similar to (2.18) is also assembled. Using (2.18) and 
the maximum spanning tree property, the branches in Au^cs and Atees are reordered 
10 in the following way. 

First, all t columns of Annks. corresponding to capacitive branches are placed 
on the left of Ah^cs- Then, from the set B^, the branches that are links to the 
capacitive trees in Guees = (iV, By) are identified and the corresponding columns put 
on the right of Auees. The final branch ordering with the corresponding 
15 permutation matrix yields a result similar to (3.7). That is 


Another TCF can be obtained by taking the RREF of (3.28) and removing zero 
rows from the bottom. This TCF has the following structure 


(3.26) 


C = size(By)-T = r-x 


(3.27) 


(3.28) 



0 hxi^"'- 0 


(3.29) 


20 


where 


T 


Z 


t 


is the number of non-capacitive link-branches in the set Bx. 
■ is the number of capacitive link-branches in the set Bx. 
is the number of capacitive tree-branches in the set By a& defined in 


(3.26). 


25 


^ - is the number of non-capacitive tree-branches in the set By as defined in 


(3.27). 


Also, ^ + T - r and t + z = s. If any of the dimensions are equal to zero, the 
corresponding block in (3.29) would simply disappear. > 
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Based on the TCF (3.29), the reduced node incidence matrix and the basic 
loop matrix, both referred to the branch ordering defined in (3.28), are found as 
Ao = Atcf (3.30) 


-fAfx»1 ; 0 0 1 


Similar to (3.11)-(3.17), KCLs and KVLs can be written for the two networks. 
That is, writing KCL using (3.30), 


(3.32) 
(3.33) 


and 


Writing KVL based on (3.31), 


(3.35) 
(3.36) 


Selecting the vector of state variables for iVc to be the vector of Independent 
voltages Vy, similar to (2.88), then 


= (Af)\ 


Expression (3.28) can also be written as 


(3.39) 
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Assuming that the vector of branch currents for network Na can be expressed in 
terms of the branch order (3.28) as 

'A 'A :a 

Ibr = GbrVbr-ibr (3.40) 

the system of linear equations for based on (3.34), (3.39), and (3.40) can be 



(of." 

1 


"A 
0 An 



If the branch conductances do not depend on voltages, the system of equations 
(3.41) has the form Ax = b, with equations and the same number of unknowns, 
and A being a full-rank square matrix. 

Using the dimensionality reduction procedure in Wasynczuk and Sudhoff, the 
state equation for the Nc in the presence of the algebraic network b^^icomes 

where the reduced quantities are defined as 
and 

Cbr = T^Ci^Ti (3.45) 


G6r = T^G,,T, 


(3.46) 


Inductive, Capacitive, and Algebraic Network Interconnection 

Heretofore, both types of network, namely A^^ and Nc, have been considered 
with a shortage of branches that can carry state variables. In both cases, the 
corresponding minimal state equations (3.18) and (3.42) were completed by adding 
extra source terms that came about due to the algebraic part of th6 corresponding 
network. Also, both types of the final systems of DAEs are stmctur^ly similar due 
to their dual circuit. These structural properties can be utilized even further in 
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obtaining a system of DAEs for the global network = (G, P) that includes Nl. 
Nc, and Na- In doing so, there are two approaches. 

In the first approach, Nl and Na are considered as was done in the beginning 
of this chapter, whereupon Nc is incorporated into the existing structure. In doing 
5 so, capacitor voltages of Nc can be mapped into e^^ and v^^ in (3.18), A second 
approach would consist of adding Nl to the structure developed for the Nc and Na, 
and mapping inductor currents into j^^and i^^ in (3.42). Since both methods yield 
equivalent results, choosing either one is a matter of pure convenience. Following 
the order in which the material was presented earlier, preference is given to the 

10 first approach. 

Suppose that the MinSTA is applied resulting in (3.5) from which the branch 
order (3.7) and the TCF (3.8) are obtained. In the most general case, Nc could 
have its branches anywhere in Gtrees = (N, By) and among any noh-^nductive link- 
branches in Giinks = (N, Bx). Relative to (3.7), these branches could be columns 

15 corresponding to blocks , A^^ , and A,'*^ . The capacitive branches 

corresponding to any of the columns in A^^ represent the part of Nc that overlaps 
with Nl- As a result of such overlapping, each of such capacitive branches is going 
to have its own independent state variable within No In fact, the branch voltages 
corresponding to such capacitive branches can be viewed as independent voltage 

20 sources e^^ present in (3.18). The remaining capacitive branches are viewed as a 
part of an algebraic network for Nl, and therefore, are going to be represented in 
the columns of blocks A^, and A^^ as given in (3.7). 

Now, the challenge is to reorder columns (branches) in (3.7) taking into 
consideration the capacitive network. The procedure of reordering columns is very 

25 similar to the two previous cases. In particular, all columns corresponding to the 
non-capacitive branches in A^^ are identified and placed on the left side of this 
block. Then, the trees of capacitive branches that form a capacitive network need 
to be separated. In order to achieve this, the MaxSTA with weight function wdb) 
is applied to the branches in A;^, and A^„^ . Then, with the result (3.26), the 

30 columns of this block are again sorted similar to (3.28). Thereafter,'''the final 
branch ordering with corresponding permutation matrix may be expressed as 


16410-112/154442 


35 


A/Tlca = [aL.A;^,aL,a1,.A^^.aSL,aL] (3.47) 
The permutation matrix Tlca in (3.47) sorts branches of the global network AT 
in groups with very specific topological properties corresponding to different 
networks. Again, taking the RREF of (3.47) and removing the zero rows, a TCF 
with the following structure is produced 


(3.48) 


i;^^ 0 0 0 0 0 A^xft 
0 I^^^ 0 0 0 0 A^x/ 
0 0 e 0 D^tA^./ci^ 
0 0 0 1^x5 Acxt 0 Cz^h 

where 

■q — is the number of non-capacitive tree-branches in the set By that are also 
placed in Nl. 

H - is the number of capacitive tree-branches in By that are placed in Nl and 

No 

T - is the number of capacitive tree-branches in By that are placed in Nc- 
^ - is the number of non-capacitive tree-branches in 5y that are placed in Na. 
f - is the number of non-capacitive link branches in Bx, that are placed in Na. 
z-is the number of capacitive link-branches in Bx. that are placed in Nc. 
h- is the number of inductive link-branches in the set Bx that are placed in Nl. 
The reduced node incidence matrix and the basic loop matrix for the global 
network referred to the branch order (3.47) are then found as 

Aa = Atcj (3.49) 

"CA r -A T 
0 0 -(Dt>cO -(Acxt) It^t 0 0 

0 0 -(A^x.)'' 0 0 I,,, 0 (3.50) 


From this point, it is possible to proceed as usual. KCL for A^z, using (3.49) 
can be written as 
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ForiVLandA^A, 

[oc><n....I^x,Afx.O,.,C^.]i*^ = » (3-53) 
[l^.,A,^]^*-[o,.,.,.,-C^:fJi^'' (3.54) 

'A'A 'LA-L 
Aalbr = -Co ibr 

Finally, KCL relating Nc to the two remaining networks is 


kx, 0,,^ I^t OtxC ktz &\ 


Based on (3.50), KVL can be written in the following way. First, jPt^r the capacitive 
network 


Then, for Nc and Na 


|o,.,..-(DS.f-{Aj«,f,,>,0,„..|v.,= « 


Finally, KVL relating all three networks can be written as 

f-{A,\.f -(I,-./ -(c-./-(cf;^.)%,.,^ = 0 (3.6.) 

[-{I^x.f 0,.,I,.J-^ - [{J:^\c^\,^l^fr (3.62) 

"^llCrxJ 0^...r'"- ^'''^tr = -Cb Vbr-Cb ftr 

As before, the vector of state variables foriVz. is selected as a vectorof independent 
currents such that > 
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Note that in (3.63), it is possible to use from (3.62) in order to avoid 
computing currents of capacitive branches. 

Similar016410-1 12to (3.38), the vector of states for iVc is chosen to be £ 
5 vector of independent capacitor voltages such that 

, = (3.64) 

Based on definition (3.63), expression (3.54) becomes 

A^i6r = -Ca'^ibr = -C"i, (3.65) ■ ■ 

Also, based on (3.63) and (3.56), it follows that 


(3.66) 


which then can be used to rewrite the compact form of (3.56) as 

~crc Z.LCL z.ca;a _tr. Ica'a 

(3.67) 


Similarly, based on (3.64), KVL (3.60) and (3.62) can be rewritten as 

Bt Vir = -Dfvtr = -Xf^Vy (3.68) 

f(^'j.)'(cf;.)'.J'£ ' \(>:f4(&i[. ■ 

BkVi,, - -C* Vfcr-C* vtr = -C-Vy-Ci, vj, (3.70) 

15 Note that KCL (3.67) written for Nc and KVL (3.70) written for Nl are coupled 

through the corresponding interconnection matrices that are related to each other as 

C^^ = -(D'Y (3.71) 

Using (3.70), (3.67), and the dimensionality reduction procedure of 
Wasynczuk and Sudhoff, the state equations for inductive and capa0itive networks 
20 become 
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(3.73) 


(3.72) 


where all reduced quantities with the subscripts "x" and "y" can be computed very 
much like the ones in (3.19), (3.20), (3.43), and (3.44) using the corresponding 
KVL matrix and KCL matrix . As before, it can be verified that (3.72) and 
(3.73) are the minimal state equations for the Nl and A^o respectively. The 
reduced quantities Lx and are indeed non-singular matrices for any practical 
network. 

As before, in order to solve (3.72) and (3.73) for the state derivatives, the 
branch currents and voltages of die algebraic network must be computed first. 
Here, it is reasonable to assume that iV^ contains non-empty resistive branches that 
are modeled as depicted in Figs. 4(a) and (b), with voltage or currept sources, 
respectively. 

Since a single branch cannot contain both types of sources and assuming that 
each branch of Na has a nonzero resistor, the voltage equation for A''^ can be 
written as 


Then, in order to obtain a system of equations for Na, KCL (3.65) and KVL 
(3.68) should be utilized. Thus, substituting (3.74) into (3.68) and combining the 
result with (3.65), the following system of equations is obtained 


which has ^ + / equations and the same number of unknowns. If R^^ does not 
depend on branch currents, the system (3.75) would have the form Ax = b with A 
being a full rank square matrix. Thereafter, (3.75) can be solved for the vector of 
currents i^^ which is then used in (3.74) in order to compute the voltages v^^ . 
Finally, (3.72)-(3.75) represent a consistent system of DAEs for the .three 

elementary networks Nl. Nc, and Na. ' •' 

'■k 

Elementary Network Identification Procedure 


V6r = R6r(iir + jbr) + e6r 


(3.74) 
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For large electrical networks, storing the associated multi-graphs in matrix 
format leads to sparse matrices of accordingly large dimensions. Even though the 
actual topological matrices used for elementary networks as KCLs and KVLs are 
smaller in size, these matrices are sparse and their dimensions are proportional to 
the size of the elementary networks. However, in order to gain computational 
efficiency, it is desirable to employ techniques that would avoid multiplication of 
large sparse matrices. A direct way to reduce the unnecessary operations is to use 
the techniques developed for sparse algebra in all computations involving 
topological matrices. Essentially this would lead to component-wise operations 
involving all of the network branch parameters. Alternatively, it is possible to use 
a more compact technique for graph representation such as collections of 
adjacency lists instead of adjacency matrices. Such lists can be conveniently 
implemented on a computer, which in turn would enable the usage of more 
efficient Min/MaxST algorithms. Representation of graphs in temiis of compact 
data structures also suggests reformulation of the network identification procedure 
using ordered sets instead of topological matrices. 
Input Data 

Given a general electrical network N, the associated graph is defined in terms of 
node and branch sets as G = (N, B). The graph G must be closed (G must form a 
valid circuit), which implies the following 

(1) each branch in 5 is a member of some cycle composed of more than one 
branch (self-loops are not allowed) 

(2) The entire branch set B can be partitioned into two subsets By and B^c such 
that they have no branches in common, set By has a minimal size and spans the 
entire node set A''. Moreover, Gtrees = (N, By) is a spanning tree. If G is a multi- 
graph, then Gt is a spanning forest (forest of spanning trees). Also, based on the 
parameter set P, the weight functions fit, (2.81) and at (2.86) must be defined over 
the entire branch set B. Then, depending on the order in which the elementary 
networks are to be found, a network identification procedure can be formulated in 
four major steps. Two procedures will now be discussed for the different order of 
network identification. ' 

L-C-A Procedure 

Step 1: Call the minimum spanning tree algorithm 
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MinSTA(G, w£) =» = (N. By) (3.80) 

Given By from (3.80), the set of link-branches is defined as 

B^=B-By (3.81) 
Branches in Bx for which oaXbj) = 1 (inductive link-branches) are identified 
defining the set . Then, the set of link-branches is re-ordered as follows 


(3.82) 


where is a temporary set of non-inductive link-branches. 
Step 2: B^ is removed from the global branch set B and the remaining closed 
graph is defined G = {n,B^'^). Then, it is necessary to find all the tree-branches in 
By that are linked by B^ . Let this temporary branch set be denoted as B^ . 
Thereafter, the branch set for G is found as 

-CA \~CA ~CA] 

B =js, ,Bx [ (3.83) 

Based on the MinST property, it can be proved that none of the branches in 
5^"* are inductive. The branch set By can be partitioned into the following sets 


Since some of the branches in B^ may be capacitive, this branch set may be 
partitioned further as 

={4^?} (3.85) 

Also note that since the branch set 5 (3.83) forms a closed graph, the remaining 
branches {5^"^,B^}do not, unless 5*^ and are sets of galvanically 

disjointed branches. 

Step 3: The maximum spanning tree algorithm is applied 

MaxSTA(G, Wc) => Gtre« = (N, By^) (^3,8$) 
Since 5^ has an optimal property, it is no longer a temporary set. .(jviso, due to 
such possible improvement, this set may be different from the one used in (3.83). 
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Proceeding further, all branches in Bf^ for which cocibj) = 1 (capacitive branches) 
are identified and placed into a separate set here denoted as By^ Then B^ ; can be 
partitioned as 

(3.S7) 

The set of link-branches is retrieved as 

B^^ = b'^-B^ (3.88) 

Step 4: By is removed from the global branch set 5j^and the remaining closed 
graph is denoted G = (n,B^). Next, all link-branches in B^ that correspond to 
tree-branches in B^^ are found. Let this branch set be denoted as 5f . This set 
actually may or may not have capacitive branches, but all of its branches are links 
to capacitive trees in B^^ . The remaining link-branches can be found as 

B^-'B^-B^ (3.89) 
Then, the following set partitioning can be formed 


5? = {<4 


(3.90) 


(3.91) 


Based on the property of MaxSTA, it is possible prove that B^ has no capacitive 
branches. 

Finally, the global branch set B can be reordered in terms of the sets 
identified so far as 

Thereafter, the elementary networks are formed based on the branch sets as 
follows " ■' ' 
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(3.93) 


The branch order established in (3.92) corresponds to the TCF (3.48) with the 
block matrices having dimensions of their respective sets in (3.93)-(3.95). 

C-L-A Procedure 
5 Step 1: The maximum spanning tree algorithm is applied 

MaxSTACG. w^) =f = {N. By) (3.96) 
The set of link-branches becomes . i 

^x-B-By (3.97) 
All capacitive branches in By are identified. This set is then re-sorted as follows 
B^ = |5j,B^^| (3.98) 

10 

where 5^ is the largest set of capacitive tree-branches, 5" is a temporary set of 
non-capacitive tree-branches. 

Step 2: All link-branches in Bx corresponding to trees in 5 " are found. This set 
of link-branches is denoted as 5 " . The combined set 

-LA f - LA -LAl 

B =iBy ,Bx [ (3.99) 

15 

has no capacitive branches, moreover (3.99) can form a closed graph 
G = (iV,5"). 

The remaining link-branches become 

= (3.100) : ' 
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Among the branches of , there may be some which are inductive. These 
inductive link-branches can be identified and placed into set denoted . After 
that, can be re-sorted as 

Bf = |B^Bf| (3.101) 

where B^ includes the remaining non-inductive link-branches that may or may not 
be capacitive. 

Step 3: The minimum spanning tree algorithm is called 

MinSTA(G, Wj) Gnees = (N, B^^) (3.102) 
The corresponding set of link-branches is 

5," = (3.103) 
In the set 5" , all inductive branches are placed in the set of inductive link 
branches denoted as B^ . The remaining non-inductive links are 

B^'-B^-B^ (3.104) 
Thereafter, the set 5 " can be reordered as 

= !^Bt,B^^ (3.105) 

Step 4: 5 ^ is removed from the branch set 5" . The remaining closed graph 
becomes G = {n,B^). In order to accomplish this, it is necessary to identify all 
tree-branches in B^ that are linked by B^ . This branch set is denoted as By . 
Thereafter, B " can be partitioned as 

< = {^^.^} (3.106) 

where B^ is found as usual 

B^^^B^-B^ (3.107) , ; 
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It is important to note that none of the branches in {b^ , }are inductive or 
capacitive. 

Finally, the global branch set B can be re-ordered in terms of the smaller sets. 
Let the order be given as follows 

B = ^B^, bJ. B^, B^, B^. B^. B^I (3.108) 

The elementary networks are formed based on the branch sets as follows 

\b^,B^A^Nc (3.109) 


|^.B^B^j6iVi 
fe.B^UiV^ 


(3.110) 


Based on the branch order established in (3.108), the TCF of the node incidence 
matrix has the following structure 


■'^TCF ' 


C CA ' In. " l-\ ^ C 

0 0 Dtxj D^xft Axxjfe Atx; 

0 ij^^ 0 0 aJxA 0 0 

0 0 l^^^M^tclth 0 0 


(3.112) 


where blocks have dimensions corresponding to the respective branch sets in 
(3.108). In particular 

T = size( By ), as determined in (3.98) 
Ti = size(By ), as determined in (3.107) 
C = size(5^ ) , as determined in (3.106) 
t = size( B^) ,as determined in (3.104) 
h = size( B^ ), as determined in (3.105) 
k = size ( Bj; ), as determined in (3.101) 
z = size( B^ ), as determined in (3.101) 
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The previous procedures are the techniques for sequential identification of 
elementary networks. Depending on the relative order in which the networks are 
identified from the global graph, the sequence of steps in the procedure may differ. 
Similar procedures may be constructed in which the order of elementary network 
identification is different from the two cases considered above. However, as it is 
expected, the results of such procedures are topologically equivalent. 
General Interconnected Networks 

In the previous chapter, electrical circuits constructed of branches with 
components selected from one of the parameter sets (2.10)-(2.13) were considered. 
The topological conditions under which a circuit of a particular kind can be 
modeled as a single network of an appropriate type were also presented. Such 
differentiation was motivated on a basis of the type of equation and state variables 
(if any) required to describe the respective network. Furthermore, it was shown 
that as more freedom is allowed in terms of the circuit components and their 
topology, it is no longer possible to describe the entire circuit as a single type of 
network. Instead, for the circuits composed of branch models Fig. 4A, it was 
necessary to partition the global network N as 

iV=JViUiV^ (3.113) 

In the case of circuits constructed from branch models Fig. 1.2 (b), the necessary 
partitioning was of the form 

N^'Nc^N^ (3.114) 

Finally, when both types of branch models Figs. 4(a) and (b) were considered at 
the same time, it was shown that the circuit can be consistently described by 
viewing the corresponding global network as being partitioned as 
JV = (iVj, u Nj^c^ Nc) u (3.1 15) 

Even though (3. 11 5) comes directly from the TCF (3.48), it may also be written in 
simplified form as 

N Nj^^j Nc^ (3.1 16) 

Several observations can also be made about the actual coupling between DEs 
(3 .72)- (3 .73). The DE (3.72) has currents as state variables, and voltages as its 
driving input. The states in DE (3.73) are voltages, and it is the currents that 
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represent the external driving force. Therefore, these two state equations are said 
to be input-output compatible. Extending this framework, models Nl and Nc can 
be viewed as interconnected dynamical systems each of which is represented by its 
own state equation that are coupled through their respective inputs and outputs and, 
5 for certain topologies, through an additional interconnection term due to the AE 
corresponding to ^Va. 
Generic network identification 

Sometimes it may also be useful to partition the global system into several 
networks based not only on branch contents but also on some other additional 
10 constraints or criteria that may come from the description of the actual system and 
the topological layout of its sections. The property that must be maintained in any 
network partitioning is the ability to automatically produce KCL and KVL 
matrices that can be used to describe the corresponding networks and the coupling 
among them. 

15 In constructing die KCL and KVL matrices, an algorithm for assembling the 

TCP is needed. Previously, in order to assemble the TCP, the multi-graph Gg = (A^, 
B) was searched for a forest of spanning trees G = {n. By ) with some particular 
topological properties reflected in an appropriate weight function. This very 
technique may be extended to accommodate other constraints in addition to those 

20 which come from branch components. For instance, it may be possible to have 
several, possibly pre-specified or guessed, networks. Then, a generic weight 
function associated with the branch set B can be defined as follows 


where is to be understood as ^-th member of the corresponding parameter set 

as defined in (2. 10)-(2. 13). If the member of is a matrix, then by convention 
Pd^} 0) is the j-th diagonal entry of the corresponding matrix. Similarly, when 
Pd^} is a vector, then the notation P^f^} (j) is the y-th element. The is the ifc-th 


for ; = L, C, A 
? = l...slze(P5) 


wlibj) = 


0. if P5U}C/) = 0 

1, ifP^UlC/)''© 
±o», if ivj 


(3.119) 


network to be identified. The ± sign in front of the infinity symbol in (3. 1 19) is 
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taken such that the branch bj € is not likely to be included 

in G*^^^ = (a^. By ) whether the search is performed for the MinST or MaxST. 


Furthermore, it is desirable to assemble separate systems of DAEs along with 
the corresponding coupling terms for some, possibly pre-specified sections of a 
larger network. Ideally, the network should be partitioned in such a way that the 
number of variables coupling the corresponding systems of DAEs is minimized. If 
such a goal is feasible and the number of state variables in each of the DEs is 
significantly larger than the number of coupling variables, the corresponding 
networks may be viewed as being weakly coupled. In terms of the topology of 
weakly coupled networks, it is reasonable to expect that the number of common or 
connecting branches is small. 

Thus, a generic step of partitioning the network N into two smaller ones 
may be considered as 


In doing so, the columns in the node incidence matrix would also be reorganized 
such that the corresponding TCP would have a particular block structure needed to 
assemble the KCL and KVL matrices on a network-by-network basis. In this 
regard, two types of structures of the right hand side of TCP have been heretofore 
encountered: lower-block triangular as in (3.8); and upper-block triangular as in 
(3.29). Even though these two TCFs were computed based on different topological 
objectives, they are in fact equivalent. Therefore, without loss of generality, for 
the network partitioning (3.120), the lower-block triangular TCP may be 
considered. Based of such structure for the TCP, the KCL and KVL matrices 
would have the following form 


(3.120) 



(3.121) 


^^ ^ 0 -{A^f I2 0 

|-(Ai)^-(A2i)=^ 0 I, 


(3.122) 


If the partitioning (3.120) is final, then (3.121)-(3.122) would already require a 
particular choice of state variables for each network. That is, if a single set of state 
variables is to be used to completely describe a network, it is necessary to have a 
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KCL or a KVL (whichever is appropriate for the selected states) to include only 
the branches corresponding to this network. With respect to (3.120), it appears 
from (3.121) that it is possible to write KCL including only the branches of iVi. 
Therefore, one can select the set of branches corresponding to the block Ai in 
(3.121) to represent states such as independent currents. Such a choice of state 
variables would require N, to be an inductive network. Relative to network A^2, it 
is possible to write an independent KVL based on (3.122). This suggests a set of 
independent voltages to be the states and the network to be capacitive. If the TCF 
with its right-hand side being block diagonal (instead of block triangular) then the 
networks corresponding to these blocks would be topologically disconnected. In 
fact, for the multi-graph = (N. B), the TCF would have as many diagonal blocks 
as there are graphs in G^. 

On the other hand, if partitioning (3.120) is not final, it is possible to 
recursively apply similar sub-network identification procedures to either or both 
iVi and iV2 each time updating the TCF. Ideally, this procedure could be continued 
until a final collection of elementary networks is obtained. The result of this 
sequential network partitioning may be symbolically expressed as 

N = nIJiNI^...kjNIkj ...^N^ (3.123) 
The corresponding TCF can be expressed in the general lower block-triangular 

I. 


A, 
A, A, 


• A„<5A„ J 


(3.124) 


^•I„,AV"' A^jfe-- 

Based on the TCF (3.124), and following the usual procedure, it is possible 
to assemble the KCL and KVL matrices for the global network as well as for each 
of the networks present in (3.123). Thus, based on (3.124), the KCL matrix for the 
k-th network (that is, the self-KCL matrix) is defined as 


(3.125) 
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Then, the corresponding KCL coupling matrix relating the k-th and r-th networks 
can be written from the k-th row of the TCF as 


[O (3.126) 


Using definitions (3.125)-(3.126), the KCL equations for all networks in (3.123) 
can be expressed as 

k k 

Khr+ Y'^aHr = 0 , for = 1, 2, ...,n (3.127) 

i= 1 

where for ^ = 1 the summation is defined to be zero. Furthermore, based on the 
structure of (3.126), it can be seen that only currents of the link branches of the 
other networks are needed in (3.127). Thus, defining new matrices as 

= ^ki (3.128) 
the KCL equations (3.127) can be rewritten as 


A^t+ SDj4t= O.for/e = 1.2, 


(3.129) 


The KVL matrices are defined in a similar way. In particular, the KVL 
self-matrix for the *-th network is determined from (3.124) to be 

The KVL matrix coupling the k-Xh. and i-th networks is assembled as 

^^'^ [-(AfA)''o] (3.131) 
Thereafter, the KVL equations relating all networks in (3.123) can be written as 


B6vt+ E B*Vt = Cforie = 1.2 rv (3.132) 

i = k+\ 


with the convention that for k = n, the summation is defined to be zero. Noting that 
only voltages corresponding to the tree branches of the other networks are needed 
in (3. 132), the following definition can be made 
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(3.133) 


Then (3.132) can be rewritten with only the tree voltages under summation as 


Bjvt+ i C^V; = 0,forfe= l,2,...,;i 


(3.134) 


i = & + i 


5 


It can also be noted that (3.128) and (3.133) are related very much like the similar 
matrices in (3.71). 

Furthermore, based on the KCL equations (3.129) and KVL equations (3.134) 
it is possible to express the state equations for general interconnected networks 

10 (3. 123). In order to do this it is necessary to assume that there is no mutual 

coupling among branches belonging to different networks. As usual, the derivation 
goes through the dimensionality reduction procedure applied to each network in 
(3.123) that has state variables. Since there are two types of such networks, 
namely inductive and capacitive, there are only two types of state equations 

15 involved. In particular, if for some k, the network type index is ^ = L> meaning that 
network Ni^, is inductive, then the corresponding state equation is 


where all the reduced matrices are defined as before in terms of . Of course, in 
20 order to solve all DEs for multiple inductive networks simultaneously and without 
additional constraints, all state equations (3.135) must be input-output compatible. 
That is, for all i in (3.135) for which the network type index is ^ = L> it is necessary 


other non-inductive networks. If it is not possible to satisfy this condition for two 
25 or more inductive networks in (3.123), then those networks should be combined 
into one. Similarly if, for some ^, ^ = C, then network iV^ is capacitive. and the 
corresponding state equation is 



(3.135) 


to have = 


0, implying that the external forcing voltage may come only from 



(3.136) 
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Here, all reduced matrices are defined in terms of A * . Also, in order to be able to 
simultaneously model multiple capacitive networks without additional constraints, 
all state equations (3.136) must be input-output compatible. That is, all external 
forcing currents should come only from other non-capacitive networks, which 
implies that for all i in (3.135) for which C = C, it is necessary to have Af = 0. 
Whenever it is not possible to satisfy this condition for some capacitive networks 
in (3.123), then those networks should be re-combined. The relation of multiple 
elementary networks is symbolically depicted in Fig. 12. There, the external 
voltage and current sources are denoted as texi and jext- Also, none of the inductive 
networks are connected with each other directly. Similarly, the capacitive 
networks are not coupled among each other through their branch variables. 
However, all networks are interconnected in such a way that each network with 
state variables has its proper inputs and outputs. On the other hand, algebraic 
networks may accommodate any inputs and outputs that are required by respective 
inductive and capacitive networks. Such an arrangement of networks allows all 
DAEs to be solved simultaneously. 

The TCF (3.124) has a very general structure applicable to interconnected 
networks. The KCLs (3.127) and KVLs (3.132) also reflect possible coupling 
among all networks in (3.123) through their branch currents and voltages. If many 
of the networks in (3.123) are mutually de-coupled or weakly coupled, it can be 
expected that many of the block matrices with double subscripts in (3.124) are 
2ero. For instance, if it is possible to make the right side of the TCF (3.124) block- 
bi-diagonal, the corresponding networks in (3.123) would be sequentially 
connected. Another useful way of formulating the relations among the networks is 
to have one (or maybe several) specific networks that represent all 
interconnections. In terms of the right side of the TCF, an attempt would be made 
to form a block diagonal structure as far down as possible by appropriately 
selecting the networks. 
Selective Network Partitioning 

When modeling a switched electrical network N=(G,P,S),it may be of practical 
interest to partition N into more than just three elementary networks Ni^ Nc, and 
Na. In general, the network partitioning may be of a more advanced form such as 
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(3.123) where the branches may be assigned to networks based on some additional 
constraints, For instance, it is reasonable to establish additional constraints such 
that some of the networks of (3.123) remain unchanged throughout the entire 
simulation study. The same goal can also be pursued on the b^is of elementary 
5 networks. 

For convenience, it is assumed that each elementary network is partitioned 
into two smaller ones as 

where the running subscript is given as ^ = I^, C, A. The network" identification in 

10 (3. 143) could be performed with the objectives given in the following subsections, 
Identincation of networks with variable parameters 

Some of the branch parameters of A^^ may depend on time as well as on 
applied currents and voltages. In this case, the equations corresponding to 
become nonlinear with time-varying coefficients. Therefore, it is desirable to 

15 partition the network as in (3. 143) such that one of the networks, say N\ , takes all 
branches with nonlinear and/or time varying parameters. It may be necessary to 
include other branches in order to make ^J. a proper network. After performing 
such networic partitioning, it may become possible to assemble two systems of 
equations for A^^ and , respectively, such that each of them posses smaller 

20 dimensions. Then, the equations corresponding to can be assembled once per 
topology, and the nonlinear equations with time-varying parameters corresponding 
toiV^ would be easier to deal with due to their reduction in size. In order to 
perform partitioning of into N \ and N ^ such that their interconnection can be 
conveniently handled using dimensionality reduction procedures, it is necessary 

25 only to assign appropriate weights to the network branches as in (3.1 19). 
Identification of switched networks 

For large switched electrical networks, the automated procedure of equation 
generation becomes quite expensive. As a result of this computational expense, 
the models of systems with fast switching, such as models that include PWM, 

30 hysteresis current inverters, converters, etc., tend to require significantly longer 
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CPU time to run the simulation. However, for a given network N = (G, S), it is 
likely that not all of its structural sections participate in the switching transients 
throughout the history of topological changes. One might, therefore, attempt to 
identify system sections with constant topology and group them into separate 
5 networks. This idea can also be applied to each elementary network of N. 

Therefore, with respect to (3.143), an elementary network iVf can be partitioned 
such that either iV^ otN^ posses a system of equations that need not be 
reassembled for each new topology. Thereafter, an attempt is made to exclude all 
such networks from the equation assembling procedures that are performed at each 

10 switching instance. Separating the global network N into its switching and non- 
switching parts would not only reduce the total amount of computations required 
per change in topology, but also provide a means for the local averaging of state 
equations for the switched subnetwork as is known in the art. 
Network Partitioning for Stability Estimation 

15 The structure of the DAEs produced by ASMG may also be utilized for the 

non-impedaiice-based stability analysis of energy conversion systems represented 
by their equivalent circuits. That is, instead of relying on the linearization of the 
state equations and using Nyquist-type criteria, it is possible to generate the DAEs 
in a form suitable for the Lyapunov analysis. For instance, if needed, a change of 

20 variables co^dh be used to rewrite the state equations in the following autonomous 
form. 

, ^^ = /^(xj4-^^(x,.Xc) (3.144) 

where giCxi,, Xc) and ^c(xl, xc) are the interconnection terms that also include the 

algebraic network. The stability of (3.144)-(3.145) can be analyzed as similar to a 
25 problem of perturbed motion and studied using Lyapunov functions. However, for 

this technique, it is necessary to establish Lyapunov functions for each system in 

(3.144)-(3.145) ignoring the interconnections to begin with. 

Alternatively, since it may be feasible to assemble a separate state equation for 

the part of the network N which is time invariant and linear, it may be 
30 advantageous to formulate the problem of stability analysis into a robust control 

framework. For instance, if it is possible to collect all nonlinear elements into a 
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separate algebraic network with no current and/or voltage sources, then the terms 
that are due to Na can be viewed as a perturbation to the rest of the system. 

The underlying power of the robust control framework is that it not only 
enables analysis of the stability of nonlinear systems in certain regions of the state 
5 space (attractors, invariant subspaces, etc.) but it also makes possible the design of 
controllers that ensure some desirable properties. In particular, considering 
quadratic Lyapunov functions, many analysis and design control problems 
associated with the model of iV can be formulated as linear matrix inequalities 
(LMIs) that can be solved numerically. In this regard, it can be very advantageous 
10 to automatically assemble state and output equations for the global network N 
directly in a form permitting a linear differential inclusion (LDI) of a particular 
type. 

REALIZATION OF THE GLOBAL NETWORK 

In this section, a global network N = (G,P) with augmented topology is 

15 considered further. As shown in the previous chapter, a network of this kind, in 
general, possesses an interconnection of the form N = Nl\uNc^Na. Such a 
network representation allows the corresponding systems of DAEs (3.72)-(3.75) to 
be assembled in a consistent manner for each elementary network with respective 
interconnections. For compactness, (3.72)-(3. 75) can be written as 

20 )) 

(4.2) 

where the global state vector x^c = [^l, Xc] consists of minimal state vectors 
25 corresponding to the inductive and capacitive networks, y ,^ = [i^, , v ] is a vector 
of branch currents and voltages corresponding to the algebraic network, and 
Ue/= [Cfcr, jfcr] is an input vector that is composed of input voltage and current 
sources. 

System (4. 1) is a first-order system of ODEs, whereas (4.2) is a system of 
30 AEs. The AE (4.2) has the form (3.23), (3.25), (3.39), (3.41) and (3.74)-(3.75). If 
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the elements of Na are some known, possibly nonlinear functions of applied 
currents and voltages, then (4.2) may not have a unique solution. In general, (4.2) 
may be solved using iterative methods for solving systems of nonlinear equations. 
Global Network with Linear Parameters 

Here, it is assumed that the network parameters are possibly time-varying, but 
do not depend on branch currents or voltages. Then, assuming that (3.75) is well- 
posed, the linear AE (4,2) always has a unique solution. Moreover, it is possible to 
solve (4.2) explicitly for the y )J . After doing so, (4.1) and (4.2) can be merged into 
one system of first order ODEs. Thereafter, a system of DAEs describing the 
global network iV can be mapped into a state space realization in its standard form 


(4.3) 


15 


(4.4) 


where the output vector is defined as y/v = [W» v^r]. 

Before, the state equations (3.72)-(3.73) can actually be written in the form of 
(4.3), some work with AEs (3.74)-(3.75) is needed. For simplicity of derivation, it 
is assumed tl^at the state variables for the Nl and Nc are chosen to be independent 
currents and. Voltages as defined in (3.63) and (3.64), respectively. Then, defining 
matrix Ma as 


BbRbr 


(4.5) 


25 which is assumed to have full rank, (3.75) can be rewritten as 


30 


The expression for the branch voltages (3.74) for can be rewritten as 


(4.6) 
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Vbr = RfcrM. 


i-Rbr 


0 


i„+R6rM. 


I + RferM . 


0 

"A 

-Bb 


ebr (4. 


a" A 
-BbRbr 


jbr = Cji^ + C^Wy + cfet + Cjjbr 


Substituting (4.6) into (3.72), and (4.7) into (3.73), and collecting terms, the 
state equations for JV^ and Nc become 

-L-(B^crcf)eJ;'-L-'cJ-cjjt 

Finally (4.8) and (4.9) can be assembled to form a state equation in its standard 
form 


■L 

LA .CA 

■c Ac 


B 


(4.10) 


The stage is set for assembling (4.4). Expressions (4.6) and (4.7) already 
provide currents and voltages for the branches of iV^. Thus, it is necessary to 
express currents and voltages for the branches of the two remaining networks, Nl 
and Nc. Starting with Nl, a vector of corresponding branch currents is determined 
from independent currents as 
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itr=(J^tfix (4.11) 
A vector of branch voltages for the Nl can be written as 

Substituting (4.11) into (4.12), an expression for branch voltages in terms of 
independent currents can be obtain as 

10 

vJ, = (R,,.£^-)(Bf)\.L,,(Bt)'5' + et (4.13) 

After substituting state equation (4.8) into (4.13) and collecting terms, the final 
15 expression for branch voltages becomes 

= t-)(B.^)^-L..(BbV(R.4'-C-C?))i. (4.14) 

-(Lj,(B^V(C''''+C^C^))v, 

20 -(V(B^)^-'(B^C^■*C^)-I^)e^^' 

where is an identity-like matrix with ones in diagonal entries corresponding to 
25 branches in Nl and zeros elsewhere. 

The equations for a, capacitive network are assembled in a similar way. 
That is, recalling (3.64), a vector of corresponding branch voltages is determined 
as 

30 v^^ = (Af)\ (4.15) 

A vector of branch currents for the Nc can be found as 
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Combining (4.15) and (4.16) yields 

5 

i?. = (G,.4f'j(Af)S.C,,(Af)'^'-J?. (4.17) 

Similar to (4.14), (4.9) can be substituted for the state derivative in (4.17). 
10 Collecting corresponding terms, the following result is obtained 

'- = (('^-V^)(^»)'-C-(A«>'c"KV^-°»"°v))v, (4.18) 
+ (C,,,(Af)''c;'(D'^'^ + DfV))i, 
15 +(Ct,(A5''c;'{Aj+D?'Df)-l'')jg' 

where l'' is also an identity-like matrix with ones in diagonal entries corresponding 
20 to branches iftWc and zeros elsewhere. 

Assuming that vectors of currents and voltages for all branches of the global 
network N can be assembled (concatenated) fix)m the corresponding currents and 
voltages of iVi, Nc, and Na as 


L , C , A 


(4.19) 
(4.20) 


30 the expressions (4. 14) and (4. 1 8) are put together to form (4.4) which becomes 
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(4.21) 


Thus, (4.10) and (4.21) form an automatically assembled minimal realization for 
the global network iV. 

10 Global Network with Nonlinear Parameters 

Heretofore, only networks with linear time-varying parameters were 
considered. Electrical networks with such parameters are linear in the sense that it 
is possible to apply superposition with respect to state variables. Such networks 
can be successfully modeled using topological algorithms and dimensionality 

15 reduction procedures and that, if needed, the resulting system of DAEs can be 
mapped into a standard form state space realization. In general, network 
parameters may not only depend on time but also be some functions of state 
variables such as inductor currents and fluxes, and capacitor voltages and charges. 
For instance, a branch resistance may depend on the current flowing through or the 

20 voltage applied to the branch. Also, if nonlinear magnetic properties are to be 
included in the model, the corresponding branch inductances become state 
dependent. This phenomena is caused by the saturation of the magnetic materials 
used in inductors. In modeling semiconductor power electronic devices, the 
effective capacitance of the junction may also depend on voltage or current 

25 depending upon the type of device considered. Therefore, it should be considered 
that 


(4.22) 
(4.23) 
(4.24) 
(4.25) 
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As before, a system of DAEs for each elementary network and their 
interconnection will be derived ignoring inactive branches for the sake of 
convenience. First, the system of AEs (3.75) for the network is rewritten as 


^^^^ 

BbRbr 


A A'' A 
-Bb -BbRbi 


jb, 


(4.26) 


which is a system of nonlinear equations that can be solved iteratively for i^, . 
After doing so, the vector of corresponding branch voltages, v^^ is computed using 
(3.74). 

To assemble DAEs for the inductive network, the voltage equation for 
branches of this network may be expressed in its most general form 


(4.27) 


20 Taking into iaccount other networks KVL (3.70) can be written as 


(4.28) 


25 


30 


where the symbol Z^^ represents all terms that are due to interconnection of Ni, 
with capacitive and algebraic networks. 

Working with nonlinear inductances, (4.24) represents an inconvenience in 
that it becomes necessary to compute the total time derivative of the matrix Lx 
which, in turn, contains partial derivatives of Lx with respect to the vector ix- A 
partial derivative of a matrix with respect to vector is a "large" and cumbersome 
object. However, the nonlinear relationship between fluxes and currents must be 
incorporated into the system of DAEs for the network Therefore, instead of 
(4.24), the relationship between fluxes and currents can be modeled using a 
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saturation function that would "penalize" the flux linkages at high currents. Using 
this technique of representing magnetic saturation, the branch fluxes are expressed 


hr = ^bAt-<?br(hr) (4-29) 


where (pbn is a vector-valued saturation function. From this point, (4.29) is 
substituted into (4.27) and the dimensionality reduction procedure using KVL 
10 (4.28) is applied to the resulting voltage equation. In doing so, the vector of 
independent fluxes is found as follows 


15 


LA-Bfcp,,[(Bj)''iJ = L,i,-(p,(iJ 


(4.30) 


where, by analogy to all other reduced quantities, (px(ix) is the reduced saturation 
function. 

A direct result of the dimensionality reduction procedure is a state equation 
with independent fluxes as state variables. For consistency, this equation is written 
20 as )) 


Tt 


RA-BteJ'.-S^^ (4.31) 


25 Equations (4.30)-(4.31) form a system of DAEs for the network iVi, Here, it is 
necessary to solve the nonlinear AE (4.30) for ix at each integration step, or more 
precisely, at each call to the derivative function (4.31). In general, iterative 
methods would be used for solving (4.30) for the vector of currents. For small 
integration steps, Newton's method with an initial guess taken as from the 

30 previous integration step could yield fast convergence. It may also be possible to 
re-parameterize the saturation function (px(ix) into {p'x(Xx). If this is achieved, then 
(4.30) would become 
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L^i^ = + (4.32) 

which is a linear equation with respect to ix. 
5 Instead of re-parameterization of the saturation function <{>xCU> ix can be 

chosen to be the state vector in order to avoid having to solve a system of nonlinear 
equations. Using the chain rule, the time derivative of fluxes in (4.30) can be 
expressed as 


dt 


'x^K^^K ^ fx _3<PxVi^ , 


dx. 


Tt 


(4.33) 


15 Combining (4.30) and (4.33), the state equation for the network Nl with currents as 
the state variables can be written in implicit form as 


K-^1^--K4^>«-bM.-X- (4.34, 


30 


The vector of branch currents is computed as usual using (4. 11), and (4.28) 
can be used to compute branch voltages. Applying the chain rule to (4.30) and 
utilizing (4. 11), the derivative of branch fluxes can be expressed as 


Combining (4.28; and (4.35), the output equation becomes 
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In (4.36), there is a term proportional to the derivative of the state vector. If (4.35) 
can be expressed in explicit form, it could be substituted into (4,36). However, 
since the derivative of the state vector is computed anyway, this term can be 
present in (4.36) without any computational disadvantage. 

A sinndlar system of DAEs can also be assembled for capacitive networks. The 
current equation for branches of Nc is 

.C pC C^qbr_jC ;•; (4.37) 


and KCL (3.67) can be written here as 


C.C ,^LA n 

-aifc^ + Z =0 (4.38) 


For the capacitive network, the nonlinear relation between capacitor charges and 
voltages can be represented in a form similar to magnetic saturation. That is 


Qfer = Cbr'^br-hA^br) (4.39) 

where <pbr is also a vector-valued function representing a part of capacitance that 
depends on the voltage. Using the dimensionality reduction procedure with (4.15) 
and (4.38), a vector of independent charges can be found as 


(4.40) 


Using the chain rule and (4.15), the derivatives of (4.39)-(4.40) with respect to 
time are found as 
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dt 


Finally, following the procedure, the state equation for Nc can also be obtained in 
implicit form as 

The vector of , branch currents is computed based on (4.3 7) and (4.4 1) as 

The equations developed thus far can be written more compactly. For the 
algebraic network, (4.26) and (3.74) have the following form 


4(it.i^Vy.etjt,i) = 0 (4.45) 


Vfer = (lA^^br, Xx,Vy, Bbn 3br, t) (4.46) 


The two state equations (4.34) and (4.43) are rewritten as 

MlCO J* = v^, vt et, t) (4.47) 


Mc(i) J* = /c(i.. ^y^iLCt) (4.48) 


where Ml(0 and Mc(r) are so-called mass matrices that can dependent on time and 
state. 

Combined, (4. 11) and (4.38) are written compactly as 
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gL{^''>K'^tt) (4.49) 


Similarly, (4.44) together with (4.15) are expressed as 


where yl= [i . v and y = [i , v are the output vectors. 

Equations (4.45)-(4.50) form a system of DAEs which describe the-global 
network iV. Here, (4.45) is a system of nonlinear equations that could be solved 
numerically. Then, (4.46) and (4.49)-(4.50) are explicit systems of AEs that are 

10 not expensive to evaluate. The implicit systems of DEs (4.47)-(4.48) may be 

solved for the respective derivatives using some efficient techniques developed for 
systems of linear equations. In such arrangements, a linear solver would be called 
at each call to the corresponding derivative function. Alternatively, there are 
efficient numerical techniques that have been designed to accommodate time 

15 dependent mass matrices. 

SWITCHED ELECTRICAL NETWORKS 

A network of a general kind was previously defined as = CG, P, s), where s 
is some given topology. If AT is a switched network, the vector s is different for 
each new topology. In order for the constructed model to be meaningful, N should 

20 be a finite eljBptrical network for each encountered topology. Throughout the 

simulation process, the network N may switch sequentially from an initial topology 
So to the final topology s^. say as 

H'^v-'Si s^,...,s^ (5.1) 

25 Based on this sequence, a so-called topology matrix can be formed as 

® = f®o.s, s;,....s. s^] (5.2) 

For each distinct topology of (5. 1), there is a minimal realization for the 
inductive and capacitive networks present in N. In general, these minimal 
30 realizations may have a different number of states for each topology. That is, for a 
given topology s,- of the sequence (5.1), the state vectors of Nl and Nc belong to 
some appropriate vector spaces 
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given topology s, of the sequence (5.1), the state vectors of Nt and Nc belong to 
some appropriate vector spaces 


* (5.3, 


(5.4) 


where the dimensions of the minimal state vectors for the networks Nl and Nc are 
denoted with respect to (5.1) as 
1^ ao.tti, a^, ...,a^,...,a^ (5.5) 

Po'Pi Pf,...,P,-,...,P^ (5.6) 

respectively. 

It is also expected that in the topological sequence (5.1), there exists an s,- and 
15 S; in which the number of state variables is maximum for the realizations of Nl and 
No, respectively. Of course, there may be several distinct topologies that require 
the same maximum number of states for their respective realizations. Assuming 
that it is always possible to obtain a non-minimal realization of a desired size by 
including some redundant states, it is possible to define the maximum numbers of 
20 state variabieS required to implement networks iVi and iVc as 

, a = max{ao, ttj, cc^, a,., a^} (5.7) 

P = max{Po.Pi....,p.,...,p.....,p„J (5.8) 

respectively. 

25 Furthermore, since for each single topology s,- the network N = (G, P, s,) must 

be finite, the same is required for the switched network N = (G, P, S). 

Symbolically this requirement can be rewritten as 

JV= (G,P,S)e A/q"^ (5.9) 

where a and are defined in (5.7) and (5.8), anda + ^q. Sometimes a 
30 and J5 are referred to as the network complexity. Also, since S is a topology 

matrix, N= (G, P. S) is equivalent to a family of networks. 
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The sequence (5.1) may contain repeated topologies, and the actual number of 
distinct topologies may be small relative to the length m of the complete sequence 
of topologies throughout the entire simulation history. In this case, a reduced 
topology matrix that includes only the distinct topological vectors may be defined. 
5 Without loss of generality, it can be assumed that in the sequence so, sj, s/, ... sj, 
... Sr, Sm for some r<m, the first r+1 vectors are distinct. Thereafter, a reduced 
topology matrix is defined as 

° (5.10) 

10 

With respect to the network JV, it is possible to say that Sr spans all encountered 
allowable topologies. 

Calculation of Initial Conditions After Switching Event 

Here, the two types of switched networks, namely inductive and capacitive are 
15 considered. For this case, it is required that the currents through inductive 

branches and the voltages across capacitive branches are continuous as the network 
undergoes topological changes. Then, the continuity conditions for inductive and 
capacitive networks can be written as 

;;; if=if..=it.and||itL<oo (5.11) 

20 ,0' vf =v^,=v^.and||v^|L<oo (5.12) 

In other words, the vectors (or more precisely trajectories) and vf^ must be 
bounded and continuous across topological boundaries. Recalling how i^^ and v^ 
are related to the vectors of independent inductor currents and capacitor voltages, 
(5.1 1)-(5.12) can also be rewritten as 

25 (Bf7iL=(BLfe'=it (5.13) 

{Af)\'^={AfJy^;=vl (5.14) 
The KVL and KCL matrices corresponding to the second topology s,^., can also be 
expressed in terms of the branch order in the respective TCF as 

Bi,=|-(ALr IL oiT^r (5.15) 
30 K,=[li, A,t, ol%"J (5.16) 
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where T^l^*' and T^^ are appropriate permutation matrices. Based on the structure 
of (5.15)-(5.16), it is possible to define the following right pseudo inverses 

(B,f;,)*=Tr[0 If., or=(BS'f (5.17) 
(Ai,)*=Tr[l£, 0 Or=(AjS'f (5.18) 
It can be noted that B^" and A^' are full-rank matrices containing only columns 
that make the basis in either case. Therefore, (5.17)-(5.18) can be used to compute 
the corresponding initial conditions for the new topology. In particular, the initial 
values of the independent currents and voltages for the topology s,.,, can be 
computed as 

v';'=Aj::^'vf (5.20) 

That is, based on branch currents and voltages just before commutation, it is 
possible to compute initial values for the state variables for the next topology such 
that the usual continuity conditions are satisfied. 

In the present embodiment, the values in the system are constrained to 
represent finite currents and voltages. Thus, in resistive networks, 

V'.i .. I^^^l- < ht i|« < ~ (5.24) 

hfL < ~, and ||vf^. , L < eo (5,25) 

In networks with state variables, this constraint may be expressed as 
. if = ifM- i6..and||4|.<- 

vf=v.';,=vf„and|vfl<. 

Non-Minimal Realization 

Thus far, all state equations for the networks Nl and Nc were assembled 
algorithmically based on topological conditions using the minimal number of state 
variables required to represent the dynamics of each network. That is, for each 


(5.26) 
(5.27) 
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individual topology, the system of DAEs constitutes a minimal state-space 
realization. Since the topology of = (G, P, s) is changing, it is reasonable to 
expect that the dimensions of the minimal state vectors may also change. In 
general, the size of the minimal realization for each network Nl or Nc may vary 
significantly throughout the sequence of topologies So, Si, s„ S;„. On the 
other hand, for analytical and/or numerical reasons, it may be convenient to have 
matrices and vectors that do not change their sizes dynamically from one topology 
to the next. For this reason, in a computer implementation of the ASMG, matrices 
Af and can be defined so that the dimensions do not change throughout the 
entire simulation ofN = (G.P,S). Whenever needed, the appropriate number of 
zero rows can be appended to the bottom in order to maintain the same matrix size. 
That is, A ^. and B ^ can be modified as 


-c 

Aa = 

A^ 


0 

~ L 

Bb = 



_0 


The total number of rows for A ^ and for B ^ is then detemunea Dy me maximum 
number of state variables (5.7)-(5.8) needed to represent networks Nl and Nc for 
any topology in the sequence so, Si, Sj, s^. Furthermore, A, f and for B [ 
using in the dimensionality reduction procedures foriViand Nc, respectively, the 
reduced matrices and Cy will be block-diagonal, with full-rank upper-left blocks 
and zeros elsewhere. Therefore, having such a structure, these matrices are also 
block-invertible. Moreover, the same block-diagonal structure is preserved under 
any non-singular coordinate transformation applied to the state variables in (2.106) 
or (2.111). 

Besides having the same state-space dimensions for all encountered 
topologies, it is also possible to make some specific state variables redundant. For 
a network AT^ with state variables and constant parameters, a minimd state-space 
realization cjgi.be assembled in its standard form 

= A^x^+B^u .:'"^5.38) 
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where <^ = L, C, meaning that the network iV; may be inductive or capacitive. 
A coordinate transformation matrix can be defined as 

(5.39) 

and the new vector of state variables as 

= K^x^ (5.40) 

where M can be any real matrix of appropriate dimensions. Note that Mx^ is an 
additional vector of redundant states. That is, any element of this vector is nothing 
more than some linear combination of the state variables in x^. Then, defining the 
right pseudo inverse of to be of the form 

K;+ = [IO] . (5.41) 

(5.38) can be rewritten in terms of the new states as 

= K^A^K^ti^ + K^B^u = A^x^ + B^u ' (5.42) 

yiv = C^K^+x^ + D^u = C^x^ + D^u 
This transformation of variables can be carried through ir, msteg^d .of (5.38), the 
more general system of DAEs (4.10), (4.2 1), and (4.45)-(4.50) is considered. 

Minimal Common Realization 

Herein, a switched electrical network iV = (G, P, S) e N^*^with the sequence 
of topologies given in (5. 1 ) is considered. It is assumed that N is switching 
between two distinct generic topologies, S/ and s,+i. As usual, for any single 
topology s, the global network = (G, P, s) possesses an interconnection of 
elementary networks of the form N ^Nlkj Nc\j Na. If inactive branches are 
included in their respective elementary networks, the commutation conditions 
(5.24)-(5.27) must also be satisfied with respect to some fixed branch order. In 
particular, for networks with state variables, the continuity of inductor currents and 
capacitor voltages across topological boundaries results in continuity of respective 
branch curreiUs and voltages as expressed in (5.1 1)-(5.12). 

In addition to the switching of the branch currents and voltages, it is also 
interesting to consider the commutation of state variables between topologies. For 
this purpose, a network with state variables iV;, where ? = L, C, is considered. If 
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is modeled using its respective minimal state equation, then it is reasonable to 
expect that the state vector will also change its dimension from one topology to 
the next. Specifically, based on (5.3)-(5.4) and with respect to Si and s,+ 1, the state 
vector X; can be related to the corresponding vector space as 

x^eSi^' (5.43) 

x^^^eSt^'*' (5.44) 

where the index y,- = Oi, fii denotes the size of the minimal state vector for 
inductive or capacitive network as given in (5.5)-(5.6) with respect to the sequence 
(5.1). Naturally, since the topologies s,- and Sf + 1, are distinct, with regard to (5.43)- 
(5.44), it is reasonable to expect that 


(5.43) 


"«*^r' (5.46) 

from which^it follows that for the simulation process of switched network iV^, the 

.1 f 

state variables are discontinuous. Note that (5.46) is likely to be the case not only 
because of. (5.45), but also due to the fact that the natural state variables 
(independent currents, fluxes, voltages, and charges) are selected algorithmically 
based only on the topological information for the active portion of the network. 
That is, the state vector x'^ is assembled based only on N = (G, P, Si). Similarly, 
the state variables in x^""' are selected based on knowledge of N = (G, P, Si+i). 
Therefore, even though the dimensions of x'^ and x^*' may happen to be the same, 
in general, the state variables will not be continuous across topological boundaries. 

As shown in the preceding section, it is always possible to obtain a non- 
minimal realization for each topology of the sequence (5. 1) such that instead of 
(5.43)-(5.44), 

x^,x^*'e9l (5.47) 
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where "y" = a, |3 denotes the maximum required size for the state vectors as 
defined in (5.7)-(5.8). Realizations of the maximum required size for which (5.47) 
holds can be obtained by appropriately applying the transformation of variables 
(5.39)-(5.41) for each new topology of the sequence (5.1). However, in general, 
5 condition (5.47) does not guarantee continuity of state variable across the 

topological boundaries, and x'^iC. 

In addition to (5.47), it is desirable to replace (5.46) with an expression of the 
form 

_ i _ i+l 

xc = x!: . (5.48) 
10 If the transformation of state variables is performed in such a way .that (5.39)- 

(5.41) holds for every switching event in the sequence (5.1), the regulting model 

would exhibit global state-space continuity, x ^ e C. The ASMG with the state 
scheduling algorithm developed in this thesis is capable of mapping each network 
incidence N = (G, P, s,) into a system of DAEs. However, there are many 

15 numerical as well as analytical reasons that make the property x ^ e C very 
desirable. 

Without condition (5.48), the simulation of the network N = (G, P, S) is 
essentially a concatenation of solutions of m+l initial value problems (IVPs) 
corresponding to the time intervals between topologies (5.1). Therefore, the 
20 simulation would requure re-computing the initial conditions and restarting the 
integration routine for each new topology. In the case of continuous state 
variables, the, need for re-initializing the integration routine would disappear. 

If conditions (5.47)-(5.48) are satisfied for all switching instances of (5.1), the 
corresponding systems of DAEs for different topologies also become compatible. 

25 That is, having state continuity x ^ e C, it is possible to apply averaging 

techniques known in the art to the state equations obtained for distinct topologies. 
Another important advantage of working with state-space realizatioiis for which 

(5.47) -(5.4&>^e enforced, is that it becomes possible to apply Lyapunov-based 
stability analysis for polytopic systems. Therefore, the development of a 

30 detenmnistic state selection algorithm subject to the continuity constraints (5.47)- 

(5.48) should be considered. 
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The first step in deriving the state selection algorithm subject to x ^ € Cis to 
define the change of variables (5.39)-(5.41) for two adjacent generic topologies s/ 
and Sj+i, and enforce (5.47)-(5.48). To achieve this, it is necessary to relate state 
vectors andx^^' through some matrices such that 

(5.49) 
(5.50) 

where ' T J^' and T J. are right and left pseudo inverses of each other, 
respectively. 

10 As done above, it is assumed that all branches, active as well as inactive, are 

assigned to their respective networks and that this particular branch order is fixed 
for both topologies. Also, for convenience of derivation, the state variables are 
r « selected to be inductor currents and capacitor voltages; whereas the dimensions of 

respective state vectors are kept constant by utilizmg (5.36)-(5.42). Then the state 

15 variables for the instance of commutation between two topologies can be related as 

1 ^ gbase .L ^ B^-CBf = - (5.51) 

= A^lvf = A]Ti(Af)\; = (5.52) 

If the appropriate redundant states are included in i , 1 '^^ and v , v 'j^' , then 
20 the transformation matrices T 'i^ and T ^ , under conditions of any proper 
commutation, can always be made to have full rank. Thereafter, 'T^' and 
T f become inverses of each other. 
The second step in deriving an algorithm for global state space continuity 
consists of applying resuhs of the first step recursively to all topologies of the 
25 sequence (5. 1). During the time interval of topology si, network N = (G, P, si) is 
modeled using a transformation of variables of the form 


i+l _ i+ Irpt i 
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where for i^ = L, C;x'^ denotes i or v ^ , and x 5 is the respective stat.e vector 
transformed by such that state variables are continuous starting from initial 
5 topology so all the way up to topology Si. This step represents an inductive 

assumption. Then, in order to keep continuity of x ^during the switching from s/to 
s,.|.i , based on (5.49), the following relation must be satisfied 

= °4'Tj^^x^''' (5.54) 

10 From (5.54), it follows that the transformation of variables for the topology s,+i is 
in fact 

where 

15 71 = 0 

Based on (5,50), and very similar to (5.54), it is clear that the inverse of ^5 can 
also be found as 

. ; ^"'k° = '^'tI'kI = jY'^'-'^Tt-'^ (5.57) 

'.' , ■' n = 0 

20 The transformation matrices (5.56)-(5.57) can be used in (5.42) in order to achieve 
X ; e C. The same transformation of variables can be applied to the more general 
system of DAEs (4.10), (4.21), and (4.45)-(4.50), in order to obtain global 
continuity of state variables. 
IMPLEMENTATION 

25 The state selection method discussed herein provides a convenient structure of 

DAEs for the ASMG. This methodology can be readily implemented on a 
computer using a matrix representation of a graph corresponding to the global 
switched network N = (G, P, S). Thereafter, all DAEs for elementary networks 
can be readily assembled based on the TCF for each topology of iV = (G, P, S). 

30 However, a brute-force implementation of the ASMG involving sparse topological 
matrices together with sparse parameter matrices would result in very long 
simulation run times in the case of electrical networks with time-varying and/or 
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non-linear parameters and a large number of branches. Therefore, an efficient 
numerical implementation of the ASMG-generated DAEs is a very important issue 
that will now be addressed. 

Instead of employing matrices for the representation of graphs m iV = (G, P, 
5 S), it is possible to use a collection of linked lists or arrays with multiple 

indirections. Running spanning tree algorithms on the resulting, more compact 
data structures requires less CPU time. Utilizing this data structure, the network 
identification procedure produces arrays containing tree and link branches for each 
elementary network. Similar arrays can also be assembled by reading the indices 

10 of the appropriate rows and colunms of the TCF matrix (3.48) or (3.1 12). 
Thereafter, it is possible to obtain arrays of branches forming a cycle 
corresponding to each link branch. Each such array represents a set of branches 
forming a loop headed by its respective link branch. Also, for each network 
branch, it is possible to assemble an array that records the loop numbers 

15 corresponding to the loops in which this branch participates. Such an array 

represents a loop participation set for the given branch. Similarly, it is possible to 
form arrays of cutset branches corresponding to each tree-branch, and arrays of 
cutset participations for all branches. These compact arrays of loop (loop 
participation) and cutsets (cutset participation) sets can be formed on a network- 

20 by-network bafeis. As it will be shown, assembling DAEs and computing their 
terms using arrays of branch sets avoids unnecessary operations and, therefore, 
significantly reduces the computational complexity. 
IMPLEMENTATION OF STATE AND OUTPUT EQUATIONS 
In order to reduce the computational effort associated with solving the 

25 network DAEs, it is instructive to recall those equations in matrix form. In this 
section, the output equations (4.13), (4.17), (4.36), and (4.44) are considered first. 
Then attention is turned to the state equations (3.72) and (3.73), and their more 
general form (4.34) and (4.43). 
Output Equations 

30 Following the sequence in which the material was presented in preceding 

sections, the inductive network is considered first. For convenience, the voltage 
equation for Ni with time- varying parameters is repeated here. 
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The voltage for the ^-th branch of Ni can be written directly from (6.1) in terms of 
summations as v,^(^) = ^ ^f.C^. Oif/O + E 01^.(0 (6.2) 

^ E L,,(^.O^it(/) + et(fe) 

5 where both indices i and k span all branches of Nl. The number of 

calculations in (6.2) can be significantly reduced if the summations are performed 
including only those branches that are coupled through mutual resistances and 
inductances. Since the coupling among the network branches is known ahead of 
time, it is possible to define sets that compactly store this information. For 

10 instance, it is convenient to define Mf as the set containing indices of all branches 
that are coupled by a mutual resistance with the k-ih branch of Nl, Also, if branch 
bk has a non-zero resistance, it is convenient to let A/f contain k as well. On the 
other hand, if b,k has zero resistance, it cannot be coupled resistively to any other 
branches and consequently Mf is empty. Similarly, Mj; is defined to contain 

15 indices of inductive branches that are coupled with bk by a mutual inductance. 

Again, Mjf contains k only if branch bk is inductive and is not coupled with other 
network branches. If bk has zero inductance, the corresponding set is 

0. Furthermore, among branches represented in there are some that 
have time-varying inductance. Thus, it is possible to identify a subset Mf c Af [ 

20 that includes all branches with time-varying inductance. Utilizing such sets, (6.2) 
can be rewritten as 

vt(fe)= X I^frr(*.0i6r(0+ Z ^"'(k.mntdn) (6.3) 

25 Assuming that sets , , and A/f are assembled for each branch of Nl, 

the vector of branch voltages v^, can be computed using (6.3) with reduced 
expense compared to (6.2). 

The same concept can be applied to the capacitive network. The output 
equation for the time-varying case of Nc is 

30 ' ■ ^ ^ C 

.C C ^dC^ c ^ '^'^^'•-if (6.4) 


16410-112/154442 


77 


Similar sets corresponding to the branches with mutual conductances A/^ and 

mutual capacitances and subsets of branches with time-varying capacitances 
Mf can be assembled for the ^-th branch of No The expression for 
computing the ^-th component of ifj. becomes 

5 S Gf,C*./)if//)* 2 J'*'(*.'«)vf,(m) (6.5) 

+ X C4,(A.a)^vf,(a)-jf,(A) 

The sets M*. M^, , and Aff, M^, Aff can be obtained for the initial 

topology and then updated for each new topology. 

This technique of summing over only relevant terms can be extended to the 
10 more general form of output equations (4.36) and (4.44). In the case of inductive 

networks, among branches represented in Mj;, there may be some branches for 

which the inductance depends on network currents. The indices of these branches 

can be placed in a set Similarly, among the capacitive branches represerited in 
, there may be some branches for which the capacitance is a function of 
15 applied voltage. These branches are collected in . Clearly, c and 

Af f cr Mf . After obtaining these additional branch sets, components of (4.36) 

and (4.44) can'^se computed as follows 


20 


25 



(6.6) 


(6.7) 


State equations 

For the purposes of numerical implementation, it is convenient to view the 
state equation, either for iVf, or foriVo in the following general form 

M(x,t)^ = A(x,i) + g(u,0 (6.8) 
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where A(x, t) represents all terms that define the state self-dynamics. The forcing 
term g(u, f) takes into account all external sources and interconnections with other 
networks. Since all inputs represented by u have the same units, the function g(u, 
0 will have the form of a summation. Finally, M(x, t) is the mass matrix. 
5 Thus, Nl with time-varying parameters is considered first. The state equation 

for this case can be written as 

The forcing vector g(u, t), with respect to (6.9), can be expressed as 

Of course, computing (6.10) by the means of full matrix multiplication is not 
efficient. Instead, it is possible to utilize the topological information about the 
matrices in (6. 10), In particular, recalling KVL matrix (3.50) and expression 
(3.62), it can be concluded that the k-th component of the forcing vector (6. 10) is 

15 the sum of external (forcing) voltages taken around the k-th loop with respect to 
the inductive network. Then, by examining the nonzero entries of the respective 
KVL matrices B , C " , and C , it is possible to assemble special sets of loop 
branches fof.e^ch term in (6.10). The set of branch indices corresponding to the k- 
th inductivd link-branch is denoted as , which contains the A:-th branch. The set 

20 £ includes only those Nl loop-branches that have non-zero external voltage 

sources. Similarly, L " contains the indices of the loop-branches that happen to be 
in the algebraic network, and L^*^ includes the loop-branches that are in No 
Clearly, L = Lj; u kj Ljf contains all loop-branches (complete basic loop) 
corresponding to the k-th inductive link-branch. These branch sets can be 

25 assembled for each new network topology. Based on these compact loop-sets, the 
*-th component of the forcing vector can be computed with reduced effort as 
follows 

/(fe) = - X Bj(/2.Z)e^,(0- X Ci''ik.m)vtim) (6.11) 
- S C^^(^,/i)v^(;i) 


16410-112/154442 


79 


The reduced inductance and reduced resistance matrices in (6.9) are defined as 
a triple product of appropriate matrices. Using to the typical matrix multiplication 
process, the entries of these matrices are found as 

L^(^'^) = '^if "^if Bj(i. m)Bia 'OLi,(m. n) (6.12) 

OT = 1 n = 1 

^.(^J) = E 2 Bfe '^)B6C;. n)Ri,im, ii) (6.13) 

m = 1 n = t 

The terms under summations in (6.12) and (6.13) are non-zero only when the 
internal indices m and n correspond to the branches that are members of the f-th 
and j-th respective loops and to the non-zero entry of the respective parameter 
matrix. Then, using loop sets earlier denoted as , each (i,j)-th element of and 
can be computed with reduced computational effort as 

■^x(^'j)= 2 2 B^(i.m)B^(/,7i)Li,(/;i,;i) (6.14) 

Rr(i.»= 2 2 B^(i.m)Bj(/.7i)Rf,(m,n) (6.15) 

For networks with time-varying inductances and resistances, (6. 14) and (6. 15) 
can be update'd by running the summation indices m and n only over those 
branches that have time-varying parameters. 

This method of computing the reduced parameter matrices is advantageous 
over performing the full triple matrix multiplication. An even better performance 
can be obtained by using loop participation sets. It can be recalled that in the case 
of Nl, only the network loop branches that are linked by inductive links make a 
contribution to (6.12)-(6.13). Thus, for each fc-th branch of Nl. it is possible to 
assemble a set that contains all network loops (loop indices) in which this branch 
participates as determined by the local KVL matrix . For convenience, the loop 
numbering can be made the same as for the rows of B ^ • It is useful to define LPk 
as the set that stores loop numbers in which the it-th branch takes part. This set can 
be assembled by recording indices of the non-zero entries of the ik-th column of 
B J . Thereafter, it is possible to express the contribution of mutual inductance 
between the m-th and n-th branches in the reduced inductance matrix as follows 


16410-112/154442 


80 


AL^(i.;) = X Z B^(i,m)Bj'C/,/OLfe,(m,;i) (6.16) 

(t e LPJ(jeLP„) 

When k = m = n, expression (6.16) defines the contribution of the self 
inductance corresponding to the ^-th network branch. In the same way, the 
5 contribution of a mutual resistance between the w-th and n-th branches in the 
reduced resistance matrix can be expressed as 

AR^(U) = 2 E B^(i,m)Bi0,n)Rlim,n) --(6.17) 

Thus, using (6.16)-(6.17), the reduced parameter matrices in the state equation 
10 (6.9) can be computed and updated in a very efficient way. 

Considering Nc with time-varying parameters, the duality with respect to the 
previous network can be utilized. The corresponding state equation is 

where the forcing vector can be expressed as 

g^(u. 0 = Afjf,-D^^it-D^^ii (6.19) 

20 In this c^se, it is possible to utilize the topological information available in the 

TCF (3.48). In particular, recalling KCL expression (3.56), it can be concluded 
that the k-th component of the forcing vector (6.19) is the sum of currents extemal 
to the capacitive network. Moreover, the non-zero entries of the KCL matrices 
a ' . D correspond to some cutset branches. The set of cutset branches 

25 corresponding to the A:-th capacitive tree branch is denoted as Cf , and the 

subset CI £ includes the branches that have nonzero extemal current sources. 
The sets Cf* and Cf^ include the indices of the cutset branches that are in the 
algebraic and inductive networks, respectively. Clearly, C,^ = C^ u Cf*u Cf^ is 
the complete set of cutset branches corresponding to the it-th capacitive tree- 

30 branch. As before, all branch sets are assembled for each new network topology. 
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Thereafter, the k-th component of the forcing vector can be computed with 
minimized effort as ^^^^^ ^ _ ^ ^^c^^ jj^c _ ^ d^C*. «)it('^) • 

JsC? , meCf'' 

5 The reduced parameter matrices can be computed using cutset participation 

sets. For convenience, the cutset numbering can be made the same as that for the 
rows of A ^ . The set CPk includes the cutset indices in which the ^-th capacitive 
branch participates. By analogy, this set can be assembled by recording indices of 
the non-zero entries of the k-th column of . Thereafter, the contribution of a 
10 mutual capacitance and/or mutual conductance between the m-th and n-th branches 
in their respective reduced matrices can be expressed as 

^C;,(i.i) = X A%m)A';U,^t)Ci,,(in,n) (6.21) 

(16 CP„)(j^ CP„) 

15 AG,ii,j)= X E A^{i,m)A%?i)G^,{m,n) (6.22) 

(ie CP„){je CP„) 

Similarly, when k = m = n, (6.21)-(6.22) define contributions of the self- 
capacitance and self-conductance of the ^-th branch of the capacitive network. 
Using (6.21)J-,(6.22), the reduced parameter matrices in (6.18) can be computed and 
20 updated efficiently. 

Building Algorithm for Networks with Variable Parameters 

After the appropriate branch sets have been assembled, the implementation of 
the inductive and capacitive networks becomes very similar. In this section, a 
pseudo code implementing the output and state equations is described with respect 
25 to the inductive network only. The same techniques and conclusions can be 
dkectly applied to the capacitive network. 

To implement output equation (6.3), the branch sets earlier defined as , 
M^, and are used. In order to have convenient and fast access to parameters, it 
is assumed that all parameter matrices use full matrix storage. If standard matrix- 
30 vector multiplication is used to implement (6.2), the computational complexity is 

(6.23) 
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However, since the density of matrices in (6. 1) is a function the number of 
mutual parameters, the complexity of (6.3) is somewhat reduced as compared to 
0 (3n^ + n). In particular, the complexity of the pseudo code in Code Block 1 
implementing (6.3) is 

which, in turn, can be anywhere between 0 (4n) and 0 (3n^ + nj, and would still 
have the same upper bound on die order of growth of 0(r^). In the previous 
expressions for complexity, n is the number of branches in the given network, 
nf,nf,n^ and nf,«^,n^ are the total numbers of self and mutual resistive, 
inductive, and time-varying inductive parameters, respectively.* The values of 
those numbers are readily determined from the branch sets M^, M\ and A/"'. 


for (VAjsjB ) { /* each A:-th branch in JV^j, * / 
/* - initialize memory slot with e^^ */ ' 
VLnw [k] = EsL {k\ ; 

/* j-th and ;c-th branches have mutual R */ 
for (SfjeM^) { 

VLnwtk] = VLnw [A:] + RLnwEic] [j] * ILnw[j]; 

/* j-th and k-th branches have mutual L */ 
for (V/e M^) { 

^ VLnw[k] = VLnwfW + Lnw[k] [j] * dILnw[j];, 

/* j-th and ic-th branches have mutual L(t) '* 
for ( y/ e iwf ' ) { 

VLnw[;c] = VLnwM + dLnwC;c] [j] * ILnwCj]; 

} 


Code Block 1 

25 Implementation of the state equation (6.8) is accomplished in parts. First, the 

forcing vector is computed according to (6.10). Second, depending on the type of 
network parameters and choice of state variables, the vector of state derivatives is 
calculated. A pseudo code block implementing (6. 10) is shown in Code Block 2. 
The code is similar to that shown in Code Block 1. The computational complexity 

30 of die pseudcTcode in Code Block 2 is 
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r^,^oin) = eU^ i X .f^+ E 4"=] ' (6.25) 

V i=i i=i i=i ^ 

which has the same upper bound of O(n^). Here, parameter n denotes the number 
5 of inductive link branches or the total number of inductive loops, and n ^ is the 
number of voltage sources in each «-th inductive loop. Then n " and are the 
numbers of algebraic and capacitive network branches, respectively, that are also 
part of i-th inductive loop. These quantities are determined by the branch sets L^, 


10 


15 


for (* =-l...size(J3^)) { 

/* each /c-th inductive link/loop */ 

gIi[W = 0.0; /* 

initialize memory slot */ 

/* j'-th branch in 

ic-th loop has voltage source */ 

for (V/e hi") { 


gL[A:] = gLCW - 

} 

BCA:] Cj] * EsL[j] ; 

/* j-ch branch in 

/c-th loop is also in iV^ */ 

for (Vye { 


gL[fc] = gL[« - 

} 

CbLA[;c] [J] * VAnw[j];' 

/* j-th branch in 

Jc-th loop is also in Nq */ 

for (V/'e L}f) { 


gL[W = gL[;!:] - 

} 

} 

C_LC[fc] [j] * vytj] ; 


Code Block 2 

20 The building algorithm for assembling reduced parameter matrices consists of 

two parts: initialization and update. The pseudo code implementing these 
procedures is given in Code Blocks 3-5. Since all reduced parameter matrices can 
be assembled using their respective loop and loop participation (cutset and cutset 
participatibiij branch sets, the building algorithm for the other reduced matrices 

25 can be obtained by appropriately modifying the pseudo code given therein. 
Therefore, only the reduced inductance matrix is discussed. 

The computational complexity as well as features of algorithms given in Code 
Blocks 3-5 are very different. The pseudo code in Code Block 3 can be used to 
assemble the reduced parameter matrices utilizing loop (cutset) sets as in (6. 14). 

30 This code ms^'be used to the initialize and update the parameter matrices. The 
initialization is performed by writing zeros in each memory slot, and then by 
computing the appropriate contribution due to the constant inductances into an 
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auxiliary storage matrix L^,. Pre-computing the contribution of all time-invariant 
inductances into a separate matrix suggests this approach 

L^(0 = L^ + AL^(f) (6.26) 


for (i = l...size(B^) ) { /* i-th inductive loop 


*/ 


for (i = l...size(S^) ) { /* j-th inductive loop */ 
LxEi] [j] = 0.0; /* initialize memory slot */ 
/* m-th branch is in i-th loop */ 

for (Vm € if) { 

/* n-th branch is in j-th loop */ 
for (VAie if) { 

Lx[i] W = Lx[i] [j] 

+ BbL[i] [m] * Bhh[j] [n] * hmt.im] [n] ; 

} 


Code Block 3 


/* inductance between m-th and n.-th branches of iV^ */" 
for ( V<7n,n,> e L^^ ) { 

/* m-th branch is in i-th inductive loop */ 
.for (VieLP^) { 

/* n-th branch is in j-th inductive loop */ 
for ( Vj G LP„ ) { 

Lx[i] ij] = Lx[i] [j] 

+ BbL[i] [m]*BbL[7] [n.] *L [m] [n] ; 

} 

} 

} 


Code Block 4 
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/* variable L between /n-th and M-th branches of iV^ */ 
for ( V<m, 7t> e L]!^^ ) { 

/* /n-th branch is in i-th inductive loop */ 
for { Vi e LP^ ) { 

/* n-th branch is in j-th inductive loop */ 
for (V/e LPJ { 

Lx[i] C/] = Lx[i] [j] + BbL[i] [m] * BbL [j] [n] 
* (Lvar [ml [n] - Lold [m] In] ) ; 

} 

} 

/* record used value for the future call */ 
Lold[m] [71] = Lvar[m] [n] ; 

} 


Code Block 5 i 

where AJUxit) is the respective contribution due to remaining variable inductances. 
Then the procedure for updating L^it) can be implemented by initializing each 
memory slot ofIjx(.t) with pre-computed values from corresponding slots in L.^. 
This technique, however, has the disadvantage of having to copy the entire matrix 
Lx at each update. In large networks where only few parameters are variable, 
copying the entire matrix L;^ over to Lx(t) may represent significant and 
unnecessary overhead. 

Utilizing the loop participation sets according to (6.16)-(6.17), the reduced 
parameter matrices can be assembled as shown in Code Block 4, without the 
disadvantage of having to copy the entire matrix for each update. Thus, the code in 
Code Block 4 avoids the scheme (6.26) altogether. However, if the code in Code 
Block 4 is used to update the reduced matrix, the original matrix will be destroyed 
which, in turn, makes it difficult to carry the contribution due to time-invariant 
parameters from one update to the next. In the algorithm shown in Code Block 3, 
this function is accomplished by copying the entire matrix; whereas Code Block 4 
updates only the relevant entries. One way of using the algorithm in Code Block 4 
for systems with both time-varying and time-invariant parameters is to "undo" the 
previous update before performing a new one. This two-step update procedure is 
equivalent to updating the matrix once with the difference between the old and the 
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new values of the variable parameters. The pseudo code illustrating the update 
procedure is given in Code Block 5. Therein, the reduced inductance matrix is 
updated due only to changes in variable inductances. The implementation of this 
update algorithm requires auxiliary static memory for storing previous update 
5 values. This storage should be of the same size and type as the original reduced 
parameter matrix. It is also noted that the subtraction of the old and new 
inductances should be performed first as indicated by the parentheses so as to 
reduce the round-off errors due to finite precision machine arithmetic. 

The complexity of the algorithms in Code Blocks 3-5 is a function of the size 

10 of the loop and loop participation sets. The sizes of these sets and their structures 
are, in turn, determined by the self parameters and the mutual coupling between 
network branches. Therefore, the complexity of algorithms Code Blocks 3-5 is 
highly systeih-dependent. Analyzing each algorithm for the worst case, which 
corresponds to all branches having variable parameters with all branches coupled 

15 with each other, results in similar complexity for all algorithms and reveals little 
about their actual performances with respect to practical networks. Instead, the 
problem of assembling the reduced parameter matrices using algorithms in Code 
Block 3 and Code Block 4 should be considered with respect to sonie typical cases. 
Expressions for the complexities can be significantly simplified using certain 

20 assumptions. Thus, for the purpose of derivation, it is assumed that there is no 

mutual coupling. Then, inspecting Code Block 3 for implementing (6.14) with this 
assumption in mind, the complexity can be expressed as 

where n is the number of inductive loops, and m is the average number of branches 
25 in each loop. Denoting m, to be the number of branches in the i-th loop, (6.27) can 
be developed further as 

n6.i4('^) = S + '^'j = ©I i: 'n.j + /I'j (6.28) 

With respect to (6.16), the complexity of the code in Code Block 4 is 
30 determined by the sizes of the loop participation sets LP. In particular, assuming 
only self inductances, the computational effort for the contribution from one 
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branch is proportional to the square of the number of loops in which this branch 
participates. Denoting the number of loops in which the i-th branch participates as 
Pi, including initialization, the expression for the complexity for the code in Code 
Block 4 becomes 


Both (6.28) and (6.29) demonstrate significant computational savings when 
compared to ®{n\nl^ +1)], which is the cost of triple matrix multiplication. Here, 
Wfcr denotes the total number of branches in the given network. A closer look at 
(6.28) and (6.29) suggests that for large networks, the advantage of the building 
algorithm presented here should become more noticeable. 

The run-time computational routines in one embodiment of the present 
invention will now be discussed in relation to Fig. 5. ODE solver (41 maintains 
state vector x, which is provided to inductive links current calculator 43 (for state 
variables relating to inductive elements) and capacitive trees voltage calculator 45 
(the portions Xc represent state variables relating to capacitive elements). 
Inductive links current calculator 43 calculates the current in the inductive link 
branches of the circuit, providing that information to resistive network algebraic 
equation calculator 47 and inductive network state/output equation evaluator 49. 
Likewise, capacitive trees voltage calculator 45 calculates the voltages v,^^^ in 
capacitive tree branches as discussed above, and provides that information to 
resistive network algebraic equation component 47 and capacitive network 
state/output (Equation evaluator 51. 

Resistive network algebraic equation evaluator 47 uses i^^^ and with e^^ 
aiid it to calculate i^, and v^, , which are provided to inductive network 
state/output equation component 49, capacitive network state/output equation 
component 51, event variable calculator 53, and the system output. Inductive 
network state^output equation component 49 uses i,^ and v^, along with inputs 

and j^^ to determine , , (provided to event variable calculator 53 and the 
system output) and dxi/dt (provided to the ODE solver 41). Capacitive network 



(6.29) 
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State/output equation component 51 uses i^^ and v^^^ along with e^^ and j^, to 
calculate i^^ , v^^ , (provided to event variable calculator 53 and the system output) 
and dxc/dt (provided to the ODE solver 41). 

The branch voltages and currents are output as vectors \br and ybr- These 
5 values are used with u^^ by event variable calculator 53 to produce event variable 
z^p, which is passed to the ODE solver 41. In addition to solving the differential 
equations for the system, ODE solver 41 monitors for negative-to-positive zero 
crossings of Zxnp- If a zero crossing is encountered, which indicates that a switch or 
switches are opening or closing, the state selection algorithm (part of state model 
10 generator 3 1 in Fig. 3) if invoked to establish a new set of state variables and to 
Update the branch sets and matrices used by the state equation building algorithm 
discussed above. 

Liductive links current calculator 43 will now be discussed in relation to 
Fig. 6. At decision block 61, it is determined whether currents or fluxes have been 

15 used to define the state variables for the system. If currents have been used, then 
the state variables from Xx, are already in the proper form and are provided as 
output. Otherwise, if fluxes have been used, it is determined at decision block 63 
whether constant or variable inductance parameters are present. If the parameters 
are constant, the output is calculated at block 67 as L~'x^ - Alternatively, if the 

20 inductance parameters are variable, block 69 determines the output vector by 
solving L^i^^ =x^. 

The anedogous calculations for capacitive tree voltage calculator 45 will now 
be discussed in relation to Fig. 7. At decision block 71, it is determined whether 
voltages or charges have been used to define the state variables for the system. If 

25 voltages have been used, then the state variables from are already in the proper 
form and are provided as output. Otherwise, if charges have been used, it is 
determined at decision block 75 whether constant or variable capacitance 
parameters are present. If the parameters are constant, the output is calculated at 
block 77 as -©^'x^ • Alternatively, if the capacitance parameters are variable, 

30 block 79 determines the output vector by solving v^^ = x^ • 
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Resistive network algebraic calculation block 47 will now be discussed in 
more detail in relation to Fig. 8. At decision block 81, it is determined whether the 
resistive network parameters are constant or variable. If they are constant, the 
output current vector is calculated as shown in block 83. If it is determined at 
5 decision block 81 that the resistive network parameters are variable, the current 
vector for the algebraic network branches is determined in block 85 by solving the 
equation shown therein. After the current vector is determined in either block 83 
or block 85, the voltage vector is determined at block 87, as shown therein. 

The inductive network state and algebraic equation determining block 49 is 

10 shown in more detail and will now be discussed in relation to Fig. 9. First, at block 
90, the inductive network branch currents i^^ and forcing term are computed. 
At decision block 91 it is determined whether the inductive network contains 
constant or variable parameters. If the parameters are constant, then the type of 
state variables being used is identified at decision block 92. For currents, the 

15 derivatives with respect to time of the inductive portion of the state vector and the 
currents in the inductive link branches are calculated. Similarly, for state variables 
defined as fluxes, those values are calculated at block 94. In either event (from 
either block 93 or block 94), the time-derivative of inductive branch currents and 
the voltages for inductive branches are determined at block 95, and the outputs are 

20 generated. 

If it is determined at decision block 91 that variable parameters are present, 
then it is determined at decision block 96 whether currents or fluxes have been 
selected as state variables for the inductive portion of the network. If currents have 
been selected, then the time-derivatives of and the inductive link-branch 

25 currents are determined at block 97. Correspondingly, if fluxes are used as state 
variables, then those values are determined at block 98. In either event (from 
either block 97 or block 98), the time-derivative of inductive branch currents, as 
well as the inductive branch voltages, are determined at block 99, and the output of 
overall block 49 is generated. 

30 Analogously to other discussion in relation to Fig. 9, the capacitive network 

state in algebraic equation component 51 will now be discussed in relation to 
Fig. 10. At block 100, values are calculated for the voltages v^^ and forcing term 
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g*^ for the capacitive network. At decision block 101, it is detemjined whether 
constant or variable parameters are present. If constant parameters are being used, 
then it is determined at decision block 102, whether voltages or charges are used as 
state variables for the capacitive network. Li the former case, the time-derivative 
5 of the state variables relating to the capacitive network and the voltages for the 

capacitive tree branches are calculated at block 103. Likewise, if charges are being 
used as state variables, those values are calculated at block 104. In either event 
(whether from block 103 or block 104), the system calculates the capacitive branch 
currents and the time-derivative of the capacitive branch voltages at block 105, 

10 and the output values are generated. 

If at decision block 101 it is determined the variable parameters are present, 
then the type of state variables (voltages or charges) is determined. at decision 
block 106. If voltages are being used, then the time-derivatives of the state 
variables relating to capacitive branches, as well as the voltages relating to 

15 capacitive tree branches, are determined at block 107. If charges are being used 
for state variables, then the equations of block 108 are solved. In either event 
(whether from block 107 or block 108), the capacitive branch currents and the 
time-derivative of the capacitive branch voltages are determined at block 109, and 
the output of overall block 51 is generated. 

20 Built-in Switching Logic 

Although, as shown in Fig. 4, the ASMG has a single switch branch, different logic 
may be specified by the user to determine when a given switch is opened or closed. 
Four built-in switch types, each implementing specific switching logic, have been 
considered in an exemplary embodiment of the invention. These switch types were 

25 selected to represent many common solid-state switching devices such as diodes, 
thyristors, transistors (MOSFET, BJT, IGBT, for example), triacs, and the like. In 
general, the built-in switching logic does not permit the opening of switches that 
would cause discontinuities of currents in inductors and/or current sources, as well 
as closing of switches that would cause discontinuities of capacitor voltages and/or 

30 voltage sourjces. Such attempts would violate Kirchhoff s current law and/or 
Kirchhoff's voltage law and, therefore, are not allowed. In some embodiments, 
violation of KCL, KVL, and/or energy conservation principles results in 
appropriate error messages. The four types of switches built into the system will 
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now be discussed in relation to Fig. 1 1, in which the switching control signal is 
denoted as Ubr. 

Unlatched Bidirectional Switch (UBS) 

The switch of this type can conduct current in either direction when the 
switching control signal Mf,r is greater than zero, and block positive or negative 
voltages otherwise. That is, w^,, > 0 implies v^, = 0 , while u^, < 0 implies 

= 0 . The switch can be opened or closed at any instant of time by controlling 
variable u^^ , subject to KCL, KVL, and energy conservation principles. 
Unlatched Unidirectional Switch (UUS) 

The switch of this type is similar to the UBS, but can conduct current only in 
the positive direction. That is, > 0 and 4, > 0 implies Uy^=6, while m^, < 0 
implies if,, = 0 . In any case, i^^ > 0 . The switch is closed when the variable 
Ci =min(Mj,,VjJ crosses zero going negative-to-positive. The sWitch becomes 
open when the variable =min(M^,,/jJ crosses zero going positive-to-negative. 
This switch should not be used to short-circuit a loop of capacitors and/or voltage 
sources, unless the resulted sum of the voltages equals to zero. 
Latched Bidirectional Switch (LBS) 

The switch of this type is also similar to the UBS with the exception that it can 
be open only at the current-zero-crossing. That is, the switch is closed when 

> 0 , which results in w^, = 0 . In this case, the event variable e, = u^^ is 
monitored for negative-to-positive zero crossing. In order to open, the switch state 
controlling variable h^^ has to become negative first. After that, the switch 
automatically opens at the next current zero crossing. This logic can be used to 
represent AC arcing switch or an ideal solid-state triac. 
Latched Unidirectional Switch (LUS) 

The switch of this type is similar to the LBS with the exception that it can 
conduct current only in positive direction. That is, u^^ > 0 and > 0 implies 
«i.r = 0 . Th^;switch is closed when the variable e, = min(u^^,v^^) crosses zero 
going negative-to-positive. The switch becomes open when the variable 

= min(M^,,/^J crosses zero going positive-to-negative. This switch can be used 
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to model an ideal thyristor. These four switch types can be advantageously 
integrated into the circuit simulation system routines. For example, switching 
analysis and topology evaluation for state selection can be optimized using the 
additional information inherent in each switch type, as will occur to those skilled in 
the art. 

Many variations on the above-disclosed exemplary embodiments will be 
apparent to those skilled in the art. For example, the processes discussed above 
can be implemented on many different hardware platforms runnmg many different 
operating systems. The computational portions might be implemented as an 
integrated portion of the overall system, or might rely on a computational platform 
such as ACSL (published by The Aegis Technologies Group, Inc., of Huntsville, 
Alabama, USA), or MATLAB/Simulink (published by The MathWorks, Inc., of 
Natick, Massachusetts, USA). 

Furthermore, in various embodiments, the processes might be'.carried out on a 
single processor or be distributed among multiple processors. 

Still further, the time-domain steps in the simulation might be the same 
throughout the system, or might vary between variables, for a single variable (over 
time), between portions of the circuit being simulated, or other divisions as would 
occur to one skilled in the art. Likewise, the numerical techniques used to perform 
integration and/or differentiation (e.g., trapezoidal, NDF, or Adam^ techniques), 
the rates and maximum error parameters can also vary and might be consistent 
across the system, vary among portions of the circuit, change over time, or 
otherwise be applied as would occur to one skilled in the art. 

Various embodiments of the present invention will also use different 
techniques to revise the state equations stored for each topology. In some cases, 
the data structure(s) that describe the state equations before a topology change 
event are modified only as much as necessary to take into account the new 
topology (i.e., only the changed portions of the circuit). Alternatively or 
additionally, new state equations may be derived in whole or in part for one or 
more topologies. In still other embodiments, a cache is maintained of state 
equations for some or all of the topologies that are encountered. A wide variety of 
caching strategies and techniques are available for use with the present invention, 
as would occur to one skilled in the art. 


16410-112/154442 


93 


The parameters of the branches in the system can also be updated during the 
simulation, using a variety of strategies. In some embodiments, a data structure 
reflecting the constant parameters is maintained. Each time the variable parameter 
values are updated in the system, the constants are copied into a new data structure, 
and the variable parameters are added to them. In other embodiments, the 
parameters for the circuit at a given time r, are stored in a data structure. When the 
parameters are being updated for processing time the variable parameters from 
time ti are subtracted from the values in the data structure, then updated for time 
ff+i, and the new values are added to the values in the data structure. 

While a detailed description of one form of nodal analysis is presented herein, 
many additional and alternative techniques, approaches, and concepts may be 
applied in the context of the present invention. Likewise, a variety of different 
circuits can be simulated using the present system, such as all-inductive circuits 
with no capacitors, circuits with resistors and voltage or current sources, but 
neither inductors nor capacitors, to name just a few examples. Parameters may be 
set for a branch or combination of branches to model ideal or real circuit 
components or sub-circuits, such as diodes, transistors, thyristors, of triacs. 

Using the techniques discussed above for time-varying parameters, multiple 
models can be used for a single physical component or sub-circuit. For example, a 
detailed, computationally intensive model might be used for a component when it 
has a rapidly varying input. Then, when the input has settled to a slower-varying 
state, a simpler, less computationally intensive model may be substituted for the 
complex one. 

Yet further, the data structures used to represent information in various 
embodiments of the present invention vary widely as well. The stated structures 
may be optimized for programming simplicity, code size, storage demands, 
computational efficiency, cross-platform transferability, or other considerations as 
would occur to one skilled in the art. 

Partitioning of branch sets as discussed herein may employ a wide variety of 
algorithms as, would occur to one skilled in the art. Some spanning tree algorithms 
that are well adapted for use with the present invention are presented in T. H. 
Gormen, C. E. Leiserson, R. L. Rivest, Introduction to Algorithms . MIT Press, 
McGraw Hill, 1993; and R. E. Tarjan, Data Structures and Network Algorithms. 
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Bell Laboratories, Murray Hill, 1983. The implementation of partitioning 
algorithms to the data structures involved will likely be a consideration for each 
particular implementation of the present invention. 

All publications, prior applications, and other documents cited herein are 
5 hereby incorporated by reference in their entirety as if each had been individually 
incorporated by reference and fully set forth. 

While the invention has been illustrated and described in detail in the drawings 
and foregoing description, the same is to be considered as illustrative and not 
restrictive in character, it being understood that only the preferred embodiments 

10 have been shown and described and that all changes and modifications that would 
occur to one skilled in the relevant art are desired to be protected. It should also be 
understood that while the use of the word "preferable," "preferably," or "preferred" 
in the description above indicates that the feature so described may be more 
desirable m some embodiments, it nonetheless may not be necessary, and 

15 embodiments lacking the same may be contemplated as within the scope of the 
invention, that scope being defined only by the claims that follow. In reading the 
claims it is intended that when words such as "a," "an," "at least one," "at least a 
portion," and the like are used, there is no intention to limit the claim to exactly 
one such item unless specifically stated to the contrary in the claim. 


